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ABSTRACT 

We develop semi-empirical models of the supermassive black hole and active galac- 
tic nucleus (AGN) populations, which incorporate the black hole growth implied by 
the observed AGN luminosity function assuming a radiative efficiency e and a dis- 
tribution of Eddington ratios A. By generalizing these continuity-equation models to 
allow a distribution P(X\Mbii, z) we are able to draw on constraints from observa- 
i-J^ ■ tionally estimated A distributions and active galaxy fractions while accounting for the 

luminosity thresholds of observational samples. We consider models with a Gaussian 
distribution of log A and Gaussians augmented with a power-law tail to low A. Within 
our framework, reproducing the high observed AGN fractions at low redshift requires 
a characteristic Eddington ratio Ac that declines at late times, and matching observed 
Eddington ratio distributions requires a P(A) that broadens at low redshift. To re- 
produce the observed increase of AGN fraction with black hole or galaxy mass, we 
also require a Ac that decreases with increasing black hole mass, reducing the AGN 
luminosity associated with the most massive black holes. Finally, achieving a good 
match to the high mass end of the local black hole mass function requires an increased 
radiative efficiency at high black hole mass. We discuss the potential impact of black 
hole mergers or a A-dependent bolometric correction, and we compute evolutionary 
predictions for black hole and galaxy specific accretion rates. Despite the flexibility 
of our framework, no one model provides a good fit to all the data we consider; it is 
particularly difficult to reconcile the relatively narrow A distributions and low duty 
cycles estimated for luminous broad-line AGN with the broader A distributions and 
higher duty cycles found in more widely selected AGN samples, which typically have 
lower luminosity thresholds. 

Key words: cosmology: theory - galaxies: active - galaxies: evolution - quasars: 
general 



1 INTRODUCTION 

Supermassive black holes reside at the center of most, if not 
all, massive galaxies. The masses of black holes are tightly 
correlated with properties of their galactic hosts, especially 
the velocity dispersions and masses of their stellar bulges 
(e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000; Haring 
& Rix 2004). If black holes grow principally by radiatively 
efficient accretion, then the statistics of quasars and active 
galactic nuclei (AGN) can be used to track the growth of 
the black hole population. One common approach to mod- 
eling this growth uses the black hole "continuity equation" 
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(Cavaliere et al. 1971; Small & Blandford 1992). Simple ver- 
sions of such models have two free parameters: the radiative 
efficiency e, which relates mass accretion rate to luminosity, 
and the Eddington ratio A = L/I/Edd, which relates lumi- 
nosity to black hole mass. With the simplifying assumption 
that all active black holes have the same fixed values of e and 
A, and that each black hole has a duty cycle (i.e., a proba- 
bility of being active at a given time) that depends on mass 
and redshift, I7(Mbh, z), one can use the observed AGN lu- 
minosity function to infer the average rate at which mass is 
being added to each mass range of the black hole population, 
thus evolving the black hole mass function forward in time. 
Model assumptions can be tested by comparing to observa- 
tional estimates of the local {z — 0) black hole mass function 
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(e.g., Salucci et al. 1998; Marconi et al. 2004; Shankar et al. 
2004; Shankar, Weinberg & Miralda-Escude 2009, hereafter 
SWM). 

The models presented in SWM made use of this simpli- 
fying assumption of fixed values of e and A. Our key find- 
ings were (i) that models with e « 0.07 and A « 0.25 yield 
a reasonable match to the local black hole mass function, 
though these specific inferred parameter values are sensitive 
to remaining uncertainties in the bolometric AGN luminos- 
ity function, and (ii) that reproducing the observed luminos- 
ity function with these models required "downsizing" evolu- 
tion, with the duty cycle for high mass black holes declining 
more rapidly with time than the duty cycle for low mass 
black holes. Additional constraints on duty cycles (apart 
from their relation to the fraction of active galaxies) can be 
inferred from AGN and quasar clustering (Haiman & Hui 
2001; Martini & Weinberg 2001), although these constraints 
depend on the assumed level of scatter between AGN lumi- 
nosity and the mass of the host dark matter halo. Reproduc- 
ing the strong observed clustering of 2 = 4 quasars (Shen et 
al. 2007) requires duty cycles close to one and minimal scat- 
ter between luminosity and halo mass (White, Martini & 
Cohn 2008; Shankar et al. 2010b). However, reconciling the 
duty cycles predicted by SWM with the lower redshift clus- 
tering {z « 1.5) measured by Shen et al. (2009) requires sig- 
nificant scatter (E « 0.5 dex) between luminosity and halo 
mass to lower the predicted clustering amplitude (Shankar, 
Weinberg & Shen 2010c). 

In this paper we consider more general continuity- 
equation models than those in SWM, incorporating a dis- 
tribution of Eddington ratios -P(A), as well as allowing for 
possible redshift evolution or black hole mass dependence of 
characteristic A values or average radiative efficiency. The 
motivation for these more general models is to make con- 
tact with a broader range of observations that offer addi- 
tional constraints on AGN and black hole evolution. Most 
obviously, many authors have now used black hole masses 
inferred from linewidth measurements or host galaxy prop- 
erties to directly estimate Eddington ratios, finding evidence 
for broad P(A) distributions that change with redshift and 
black hole mass (see Section [3 . 31 below for observational ref- 
erences). Active galaxy fractions provide another set of im- 
portant constraints on models. In single-A models the duty 
cycle C/(Mbh, z) — the fraction of black holes of mass Mbh 
at redshift z that are active at a given time — follows sim- 
ply from the ratio of the observed luminosity function to 
the calculated black hole mass function (see SWM). With 
a broad P(A), the definition of duty cycle depends on the 
adopted threshold in luminosity or Eddington ratio, and be- 
cause observational samples inevitably include these thresh- 
olds, the observed active galaxy fraction also tests models of 
P(A). Other authors have recently implemented continuity 
equation models with broad Eddington ratio distributions 
(Merloni & Heinz 2008; Yu & Lu 2008; Cao 2010), but their 
studies are limited to a few specific models. Our implemen- 
tation here draws on some of the techniques introduced by 
Cao (2010) and also on techniques and ideas from Steed 
& Weinberg (2003), who attempted a broad-ranging study 
of black hole evolution within a general continuity equa- 
tion framework. The broad -P(A) models described here also 
provide a more general framework for modeling quasar and 
AGN clustering, a topic we will examine in future work. 



Phenomenological models of the black hole population 
like those in this paper complement models based on hydro- 
dynamic simulations (e.g., di Matteo et al. 2005) or semi- 
analytic models that tie quasar activity to major mergers of 
galaxies or dark matter halos (e.g., Kauffmann & Haehnelt 
2000; Wyithe & Loeb 2003; Hopkins et al. 2006). The nu- 
merical and semi-analytic approaches adopt a more com- 
plete physical scenario for the mechanisms that drive black 
hole accretion, and they tie AGN to the underlying galaxy 
population by construction. However, it can be difficult to 
interpret the significance of discrepancies with observations, 
or to know whether observational successes imply that the 
physical assumptions are correct. Phenomenological models 
are free to be driven by the data, within the constraints 
imposed by the adopted parameterization. 

The next section describes our methods for computing 
black hole evolution and duty cycles, and it summarizes the 
SWM estimate of the bolometric luminosity function that 
is our basic observational input. All of our models repro- 
duce this luminosity function by construction. In Section [3] 
we describe our three basic models for P(A) — a 5-function, 
a Gaussian in log A, and a Gaussian plus a power-law tail 
to low A — then examine their impact on mass function 
evolution and conduct a first comparison to observed A dis- 
tributions. In Section [4] we turn our focus to active galaxy 
fractions and discuss how data might favour redshift evo- 
lution and possibly mass dependence of P(A). In Section [S] 
we revisit the comparison to observed A distributions, the 
impact of black hole mergers at the rate predicted by hier- 
archical galaxy formation models, the relation between spe- 
cific black hole growth and specific star formation rate, and 
the impact of adopting a A-dependent bolometric correction. 
In Section [S] we discuss our results highlighting some inter- 
nal tensions between data sets and key tensions between 
models and data. We conclude in Section [T] The appendices 
present some details of our calculational methods and of 
our estimates of observed active galaxy fractions. Through- 
out the paper we use cosmological parameters flm ~ 0.30, 
r^A = 0.70, and h = //o/lOO km s"^ Mpc"^ = 0.7, and we 
compute the linear matter power spectra running the Smith 
et al. (2003) code with F = 0.19, and as = 0.8. 



2 FORMALISM 

2.1 Evolving the black hole mass function 

To study the global evolution of the black hole population 
through time, we develop models in which the black hole 
mass function is self-consistently evolved via the continuity 
equation. 



dt 



-(MBH,t) 



a((A/BH)nBH(MBH,t)) 



(1) 



where (Mbh) is the mean accretion rate (averaged over 
the active and inactive populations) of the black holes of 
mass Mbh at time t (see, e.g., Cavaliere et al. 1971; Small 
& Blandford 1992; Yu & Tremaine 2002; Steed & Wein- 
berg 2003; Marconi et al. 2004; Merloni 2004; Yu & Lu 
2004; Tamura et al. 2006; SWM; Shankar 2009; Raimundo 
& Fabian 2009; Kisaka & Kojima 2010). This evolution is 
equivalent to the case in which every black hole grows con- 
stantly at the mean accretion rate (A^bh). In practice, indi- 
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vidual black holes turn on and off, but the mass function evo- 
lution depends only on the mean accretion rate as a function 
of mass. Note that Eq. ^ neglects any contribution from 
black hole mergers, which do not add mass to the black hole 
population but can alter the mass function by redistribu- 
tion. In Sect ion [5. 3 1 we discuss the impact of including black 
hole mergers at the rate predicted by hierarchical models of 
structure formation. 

We first define the Eddington ratio, 



A : 



where 



LEdd = IMbh = 1.26 X 10'' 



Mb 



Me 



erg s 



(2) 



(3) 



is the Eddington (1922) luminosity computed for Thom- 
son scattering opacity and pure hydrogen composition. The 
growth rate of an active black hole of mass Mbh with Ed- 
dington ratio A is Mbh —AdBu/tci, where the e-folding time 
is (Salpeter 1964) 



4 X 10" 



yr ■ 



with 



/ — e/(l — e) and e — 



Minflo 



(4) 



(5) 



The radiative efficiency e is conventionally defined with re- 
spect to the large scale mass infiow rate, but the black hole 
growth rate Mbh is smaller by a factor of (1 — e) because of 
radiative losses. 

We then define the probability for a black hole of a 
given mass Mbh to accrete at the Eddington ratio A per 
unit log A, at redshift z as P{X\A4bii, z), normalized to urnty, 
i.e., / dlogAP(A|MBH,z) = iQ This is the quantity that 
can be related to the observed Eddington ratio distributions 
discussed in Section [T] The average growth rate of all black 
holes can be computed by convolving the input P(A|AfBH, z) 
with the duty cycle [/(Mbh, z), 



{Mbh)= / dlogAP(A|MBH,^)A[/(A/BH,^; 



Mb 



(6) 



where the integral extends over all allowed values of A. A 
physically consistent model must have U{X, Mbh, 2) ^ 1 for 
all Mbh and z. In models where the Eddington ratio has a 
single value Ai, P(X\Mbh, z) = Ai(5(A — Ai), the duty cycle 
is simply the ratio of the luminosity and black hole mass 
functions0 



f/(MBH,2) = 



HL,z) 



<E>bh(Mbh, z) 



L = A/Mbh . (7) 



^ Throughout the paper, log denotes base-10 logarithm and In 
denotes natural logarithm. 

^ Following SWM, we will always use the symbol $(x) to denote 
mass and luminosity functions in logarithmic units of L or Mbhi 
i.e. $(a:)dlogx = n{x)xln{10)dx, where n{x)dx is the comoving 
space density of black holes in the mass or luminosity range x — > 
X -\- dx, in units of Mpc~^ for h = 0.7. Unless otherwise stated, 
all masses are in units of Mq and all luminosities are bolometric 
and in units of ergs~^; e.g., logL > 45 refers to quasars with 
bolometric luminosity L > lO'^^ergs"^. 



This paper presents predictions from several models 
that take as input the observed luminosity function ^(I/, z), 
an assumed Eddington ratio distribution of active black 
holes, P(A|Mbh, z), and an assumed radiative efficiency e, to 
solve Eq. ^ forwards in time. At each timestep, we use the 
black hole mass function $bh(Mbh, 2) derived from the pre- 
vious timestep (or from an assumed initial condition in the 
first timestep), and we compute the duty cycle [/(Mbh, 2) 
by numerically solving the equation 



HL,z) 



dlog AP(A|Mbh,z)[/(Mbh,2)'1>bh(Mbh, 



(8) 

with Mbh = L/{\1). We then compute the average growth 
rate as a function of black hole mass via Eq. and, finally, 
update the black hole masses and black hole mass function 
via the continuity equation (Eq. [1} . 

Details of our numerical implementation are given in 
Appendix|X] Directly solving Eq. (|8} for U{Mbu, z) is feasi- 
ble for a discrete distribution of Eddington ratios. However, 
for all of the models discussed in this work, we will adopt 
continuous P(A|Mbh,2) distributions, which can be more 
efficiently handled following the method described by Steed 
& Weinberg (2003) and Cao (2010). The latter is based on 
solving Eq. ((8)1 for the full mass function of active black 
holes 



A^act(MBH,2) = $Bh(Mbh,2) X [/(Mbh,^) ■ 



(9) 



More specifically, at any redshift z a parameterized ac- 
tive mass function A'^act(MBH, 2) is derived by directly fit- 
ting the AGN luminosity function ^(L, z) at the same 
redshift via Eq. ((Sjl. The function A'^act(A'/BH, 2) is then 
inserted in Eq. ([T} by replacing (AfBH)'^BH(MBH, i) ~ 
/dlogAP(A|MBH,2)MBHA^act(MBH,2)/ts, thus allowing 
the computation of the full black hole mass func- 
tion $bh(A^bh,2) and the duty cycle [/(Mbh, 2) = 
iVact (A/bh , 2)/-I'bh (Mbh , 2) . 

Eq. lU is solved by imposing an initial condition that 
we fix, as in SWM, by assuming a constant value of the ini- 
tial duty cycle. We usually choose this to be 0.5 at 2 = 6 for 
all masses, but in some models we lower it to 0.1 to allow 
the duty cycles at later times to be always lower than unity. 
The results at later redshifts are insensitive to the initial 
conditions for mass bins that have experienced substantial 
accretion growth (Marconi et al. 2004; SWM). Note that 
Eq. IT]) could be solved backwards in time as well by forcing 
the initial condition to be the local black hole mass func- 
tion (e.g., Merloni 2004). In this case, incorrect assumptions 
about the physical parameters or present day mass function 
manifest themselves as unphysical mass functions at earlier 
redshifts (e.g., negative space densities). Given the still sub- 
stantial uncertainties in the local mass function (e.g., SWM; 
Vika et al. 2009; Shankar et al., in prep.) and our desire to 
explore a wide range of physical models, we prefer to evolve 
forward in time and compare to the local mass function as 
one of several observational constraints. 

^ As discussed in Appendix|Xl we adopt a double power-law form 
for the input Afact(MBH, 2), which is a good (though not perfect) 
approximation of the functional form adopted for the AGN lu- 
minosity function (see Section 12.21 1. Our procedure allows us to 
reproduce the AGN luminosity function at the ~ 5% level, ade- 
quate precision given the statistical uncertainties in the data. 
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Figure 1. Left: Comparison among the three P(A) distributions taken as a reference throughout the paper. The dotted line is the 
(5-model centered on logAc= —0.6, the solid line is the G-model, a Gaussian with dispersion of S;^= 0.3 dex and centered around the 
same value of log Ac= —0.6, and the long-dashed line is the G+P-model, characterized by the same Gaussian plus a power-law with slope 
a = —0.3 normalized to have the same value as the Gaussian at log A = log Ac — logS^j. Note that we define our distributions in log A 
but refer to them in figures and text as P(A), and that all distributions are truncated at A = 1. Right: Dimensionless emissivity per log A 
interval, AP(A), for the G and G+P distributions. Note that even the G+P model has its emissivity peak near the peak of the Gaussian 
component. 



2.2 Input AGN Luminosity Function 

A full discussion of our adopted AGN luminosity function 
appears in SWM; here we provide a brief summaryQ We 
adopt the Ueda et al. (2003) fit to the number of sources per 
unit volume per dex of luminosity logLx = logL2-iokcV, 
composed of a double power-law and an evolution term: 



^{Lx,z) = e(z,Lx) 



A 



(10) 



where 



e{Lx,z) = 
with 

Zc{Lx) ■- 



[l + zY^ 



1 + z 



\i z < Zc{Lx) 
if 2 ^ Zc{Lx) 



z',{Lx/La)° 



if Lx > La 
if Lx < La 



(11) 



(12) 



The values of the parameters have been tuned to correctly 
match the full set of data presented by SWM (see their Ta- 
ble 1). 

The X-ray luminosity function of Eq. (|10p is con- 
verted into a bolometric luminosity function using the fit 
to the bolometric correction given by Marconi et al. (2004) , 
log L/Lx = 1.54 + 0.24C + 0.012(2 _ o.o015C^ with C = 
log L/Lq — 12, and Lq — 4x 10^"^ ergs~^. We have changed 
our procedure slightly with respect to SWM, dropping the 
convolution with 0.2-dex Gaussian scatter in bolometric cor- 
rection, because the assumption that this scatter is uncor- 
related with X-ray luminosity seems insecure. However, this 
change makes minimal difference to the derived luminosity 
function, and our results using the original SWM fit would 
not be noticeably different. 



* The full tables of the input AGN luminosity functions, and 
also black hole mass functions and duty cycles (for a set of 
representative models), can be found in electronic format at 
|http : //mygepl ■ obspm. fr/~f shankar/ 1 



The luminosity function data that we fit also include 
a number density of Compton-thick sources in each lumi- 
nosity bin equal to the number of sources in the column 
density range 23 ^ log Nh /cm~^ ^ 24, computed following 
the Ueda et al. (2003) prescriptions for the P{Nh\L) column 
density distributions. We checked that our estimates of the 
Compton-thick number densities at high redshifts are fully 
consistent with the recent results by Alexander et al. (2011) 
and Fiore et al. (2011). As extensively discussed by SWM, 
the integrated intensity obtained from our model AGN lumi- 
nosity function is also consistent with all the available data 
on the cosmic X-ray background for energies above 1 keV. 

Our qualitative conclusions would not change if we 
adopted the Hopkins et al. (2007) luminosity function in 
place of SWM, though best-fit parameter values would be 
somewhat different. 



3 EFFECTS OF A BROAD P(A) 

In this section and the ones that follow, we will compare 
predictions of our models to a variety of observational con- 
straints. If the observational uncertainties were described 
by well defined statistical errors, then the natural approach 
would be to determine best-fit model parameters via min- 
imization and compare different models via Ax^- We have 
taken this approach in some of our previous studies when 
we were focusing on a specific observational constraint and 
a limited class of models. However, the uncertainties in the 
observational constraints we consider — AGN luminosity 
functions, the local black hole mass function, distributions 
of Eddington ratios, and duty cycles estimated from active 
galaxy fractions — are in all cases dominated by system- 
atic errors, and in some cases even the rough magnitude of 
these systematics is difficult to estimate. We discuss these is- 
sues in the text and Appendices below, and the papers cited 
for each observable often discuss systematic uncertainties at 
length. 

Given this situation, we adopt a philosophy of "qual- 
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itative comparison to quantitative data". We infer model 
parameters based on an overall match to one or more sets 
of observational constraints, but we do not attempt formal 
X^-minimization because the errors themselves are not well 
defined enough to do so. We note where models fit obser- 
vations within a plausible range of systematic uncertainties, 
and where they do not. Our objective is to delineate the 
global characteristics that successful black hole accretion 
models must possess to reproduce the range of observables 
now emerging from deep, large area surveys of the AGN and 
galaxy populations. 

3.1 Input Eddington ratio distributions 

In SWM we discussed accretion models for black holes with 
a single value of A, examining the effects of redshift and mass 
dependence of this value of A and describing the impact of 
uncertainties in the input AGN luminosity function related 
to its observational determination, obscured fractions, and 
bolometric corrections. From the match to the local black 
hole mass function, we were able to set constraints on the 
average input radiative efficiency and evolution in the char- 
acteristic A. In this section we study the impact of broaden- 
ing the input P(A) distribution. To this purpose we adopt 
and compare the outputs from three different distributions 
(shown in Fig. [ij : 

• the 5-model, centered on log Ac= —0.6 (dot-dashed line 
in panel a; we choose this particular value because it pro- 
vides a good match to the local black hole mass function); 

• the G-model, a Gaussian in log A with dispersion of 
Ea= 0.3 dex and centered around the same value of log Xc— 
—0.6 (long-dashed line in panel a); 

• the G+P-model, characterized by the same Gaussian 
plus a power-law with slope a = —0.3 normalized to have 
the same value as the Gaussian at log A — log Ac — log Ea 
(panel b). 

The G-l-P distribution has a shape close to the one inferred 
by Kauffmann & Heckman (2009; see also Aird et al. 2011) 
in the local Universe from the Sloan Digital Sky Survey 
(SDSS), and it is also consistent with some theoretical and 
semi-empirical expectations (e.g., Merloni & Heinz 2008; Yu 
& Lu 2008; Hopkins & Hernquist 2009; Shen 2009). In par- 
ticular, the power-law component could represent the effect 
of either a steady decline in the accretion rate after a near- 
Eddington growth phase or a second mode of AGN fuel- 
ing triggered by secular instabilities instead of major merg- 
ers. As discussed by SWM and demonstrated further below, 
the chosen value of log Ac produces reasonable agreement 
between the local and accreted mass function, especially 
around (1 — 3) x 10* A!q where the former is best determined. 
Our chosen value of a, on the other hand, was mainly de- 
rived a posteriori after some trial and error. Higher values of 
a can induce unphysical duty cycles U{Mbu, 2) > 1 in some 
mass bins during the evolution, while lower values give re- 
sults that are not much different from the G-model alone. 
For all models we adopt a radiative efficiency of e = 0.06, 
unless otherwise stated, within the range of values favoured 
by SWM. In future sections and figures we will consider 
many variants on these three basic models. For reference, 
we provide a full list of models and their basic properties 
in Table [1] We summarize the comparison between our six 



primary models — the three introduced here, plus versions 
that introduce redshift and mass dependence of P(A) — and 
observational data in Table [5] (Section 

The left panels of Fig. [2] illustrate the effect of vary- 
ing the input Eddington ratio distribution P(A) at fixed 
duty cycle U (Mbh , z) and fixed black hole mass function 
'l>BH(AfBH, z), at z = 2 and z — 1. The luminosity function 
in each case is computed via the convolution of the respec- 
tive P(A) distributions with the product of the duty cycle 
and black hole mass function (Eq. |8]). In the case of the 
5-model, the match to the input AGN luminosity function 
(Section l2.2|l is exact. For these panels we fix U{Mbu, z) and 
^bu{Mbh, z) to those predicted by the 5-model (evolved in 
accord with the observed luminosity function), so that we 
can isolate the impact of varying the input P(A) distribu- 
tion. Switching from a 5-function P(A) to a Gaussian of 
Ea = 0.3 dex has a mild impact on the predicted luminos- 
ity function. However, adopting the G-l-P distribution at 
fixed (7 (Mbh, 2) drastically lowers the luminosity function 
at high L, since the probability is dominated by low Edding- 
ton ratios. The shape of ^{L, z) parallels that of the other 
models above the break luminosity, but the G-l-P model has 
a steeper faint-end slope reflecting the contribution of low-A 
activity. 

The procedure we follow in the rest of this paper is to 
compute U{Mb¥i, z) for a given input P(A) distribution in 
a way to reproduce the observed AGN luminosity function. 
The right panels of Fig. [2] show, at the same redshifts 2 = 2 
and 2 = 1, the duty cycles inferred by matching the SWM 
luminosity function when adopting the three -P(A) distri- 
butions and the same underlying "l>BH(AfBH, 2), again that 
of the 5-model. The duty cycles I7(Mbh,2) for the 5 and 
G models are similar, but they are several times larger for 
the G-l-P model at high masses. As described in SWM, the 
duty cycles must decrease with mass in order to reproduce 
the "downsizing" effect in AGN evolution (the characteris- 
tic AGN luminosity increases with redshift). In the G-l-P 
model the duty cycle can reach values close or around unity 
at 2 < 2 and masses ~ IO'^Mq, although with low values 
of A for the majority of them. The decline of the duty cy- 
cle in this model at masses lower than Mbh ^ Lt/{l\c), 
where L* is the break of the luminosity function (Eq. |10)). 
is induced by the presence of many AGNs radiating at low 
A from high mass black holes, which can already account 
for the low-luminosity AGNs. However, the detailed form 
of this decline is sensitive to the precise shape of the AGN 
luminosity function and our assumed -P(A). 

3.2 Global Accretion Histories 

Fig. 13] shows the evolution of the black hole mass functions 
and duty cycles for the 5, G, and G-l-P models, respectively. 
In the upper panels, grey bands show the estimate of the 
local black hole mass function by SWM. The width of this 
band already encompasses a number of systematic uncer- 
tainties, but, as SWM discuss, the inferred mass function 
for Mbh < 10* M© is sensitive to uncertainties in the treat- 
ment of spiral galaxy bulges. Vika et al. (2009) have tried to 
address this problem by estimating the local black hole mass 
function on an object by object basis, i.e., by computing the 
bulge fraction for each galaxy in a large sample, assigning 
black hole masses from an Mbh-L relation, and then com- 
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Figure 2. Left panels: Predicted luminosity function at 2 = 2 (top) and at 2 = 1 (bottom) for the tliree P(A) distributions of Fig.[T] as 
labelled, when the same black hole mass function '3?bh(^bHi2) and duty cycle U{Mbii, z) are used, in this case the one predicted by 
the 5-model. Points show the data compilation of SWM. Right panels: predicted duty cycle (7(Mbhi z) at the same redshifts 2 = 2 and 
2 = 1 when, instead, the AGN luminosity function is fixed to the one in SWM, as in our self-consistent models. Allowing a large fraction 
of sub-Eddington accretion rates, as in the G+P model, increases the probability and thus the duty cycle for more massive black holes 
to be active at different luminosities. 



Reference Models 


shape & properties 


Section /Equations 


5 


5- function 


Sec.l3l1 


G 


Gaussian 


Sec. 1311 


G+P 


Gaussian+Power-Law 


Sec.lXn 


G(z) 


Gaussian+Ac(2) 


Eq.[l5] 


G(z,AfBH) 


Gaussian+Ac(2, Mbh)+1ow-2 broadening+mass-dependent e 


Eas.[T5l[T6l[T7l[T9l 


G+P(z,A/bh) 


Gaussian+Ac(2, AfBH)+2-depondent Power-Law+mass-dependent e 


Eas.[T5l[T6l[T8l[T9l 


Test Models 


shape & properties 


Ref. Eqs. 


G(z,MBH)+constant e 






G(z,MBH)+e(z) 


G(z,MBH)+z-dependent e 


Eq.[20l 


G(z,MBH)+i^(A) 


G(z,A/BH)+A-dependent bolometric correction 


Ea.[22l 


G(z,Mbh)+P(2 < 0.7) 


G(z,A/BH)+Power-Law only at 2 < 0.7 





Table 1. List of models explored in this work along with basic explanation and Sections and Equations in which they are first introduced. 
The models are divided into two groups, the "Reference Models", i.e., the three ones introduced and discussed in Section l3.ll and the 
"Test Models" which are simply variants of the Reference Models and progressively introduced in the rest of the paper. The Gaussian 
always has dispersion of logS^ = 0.3 dex except for the models with broadening, for which logS^ steadily increase with decreasing 
redshift as in Eq. II17I I. The redshift and mass dependence in the characteristic Eddington ratio Ac are defined in Eqs. JTSj and III6I I. If 
no evolution is assumed then log Ac = —0.6. The power-law component of the Eddington distribution is always normalized to have the 
same value as the Gaussian at log A = log Ac — log and a slope in log A of a = —0.3. In the last quoted model the power-law has slope 
a = —0.9 and is normalized to have the same value as the Gaussian at log A = log Ac — 0.21ogS;^. The input radiative efficiency can be 
constant or vary with mass or redshift (see Eqs. I19l and l20l l as detailed for each model. 
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Figure 3. The left, middle, and right columns refer to outputs of the 5, G, and G+P models, respectively, whose P{X) distributions are 
shown in Fig.[T] Upper panels: Predicted black hole mass function at 2 ~ {solid lines) compared to the local estimates by SWM (grey 
area), and Vika et al. (2009; squares with error bars, with the open symbols indicating the estimates derived from galaxies below their 
reliability limit). Other lines show the predicted black hole mass function at earlier redshifts, as labelled. Bottom panels: Corresponding 
duty cycles (for all sources accreting above any A > 0) as a function of black hole mass at different redshifts, as labelled. For all models 
a radiative efficiency of e = 0.06 has been assumed. 



puting the black hole mass function applying the X^/VInax 
method. Their result, shown with open and filled squares in 
Fig. O agrees well with SWM at high masses, while it turns 
over at M < 10* Mq rather than continuing to rise to lower 
masses. We thus consider Vika et al.'s and SWM's results 
to broadly bracket the still remaining uncertainties in the 
determination of the local black hole mass function, and we 
will use both as a reference when comparing with models. 

Moving from the 5-model to the G-model makes almost 
no difference to the evolution of the black hole mass function 
or to its 2: = value. Both predictions are insensitive to our 
assumed initial conditions at 2 = 6 because even by z = 4 
the accumulated mass from accretion greatly exceeds that 
in the initial seed population. The G+P model differs in this 
regard because reproducing the z = 6 luminosity function 
even with a duty cycle of 0.1 requires a high space density 
of seed black holes, since many black holes are active at low 
A values that do not contribute to the observed range of 
the luminosity function. Once the evolved black hole mass 
function substantially exceeds the seed population, evolution 
is similar to that of the 6 and G models. 

The clustering of quasars is a diagnostic for the space 
density of the underlying black hole population, as more 
numerous black holes must reside in lower mass halos that 
are less strongly clustered (Ifaiman & Hui 2001; Martini & 
Weinberg 2001). Applying this idea to the strong clustering 
measured for z f» 0.4 quasars by Shen et al. (2007) in the 
SDSS, White et al. (2008) and Shankar et al. (2010b) find 
that duty cycles close to unity are required assuming that 
these highly luminous quasars have A > 0.1. This finding 



is at odds with the high black hole density required for the 
G+P model at high redshift. To quantify this statement, we 
have computed the large scale bias of the G+P model using 
a formalism that we will describe more fully in a future pa- 
per. In brief, we match black holes to halos and subhalos via 
cumulative number matching (similar to Conroy & White 
2012), including 0.3-dex scatter in black hole mass about 
the mean relation (see Appendix |B|, then assign each black 
hole a luminosity based on the duty cycle and the Edding- 
ton ratio distribution P(A|A/bh,2). We can then compute 
the mean large scale halo bias for black holes in a luminos- 
ity range. For the G+P model, which has high black hole 
space density and large scatter between luminosity and halo 
mass, we predict a bias of 7.5 for quasars with logL > 47 at 
z — 4, which is more than 2a discrepant with the reference 
value of 12.96 ± 2.09 given by Shen et al. (2009), though it 
is marginally consistent with their lowest estimate. We con- 
clude that broad P(A) distributions at very high redshifts 
are observationally disfavoured, a conclusion in qualitative 
agreement with the findings of Cao (2010). If -P(A) does have 
a long tail to low A values, it must develop after the earliest 
stages of black hole growth; we will consider models with 
growing power-law tails in Section [4] 

Duty cycles U (Mbh , z) tend to decline with redshift 
and with mass, for the reasons already explained in Fig. 2. 
We note again the higher values of the duty cycle in the 
G+P model at all redshifts. The high-redshift curves are still 
affected by the initial conditions in the black hole population 
at low masses. 

Fig. [4]tracks the overall growth of the black hole popu- 
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Figure 4. Upper panels: Growth rate of the integrated black hole mass density as a function of time {solid lines) and relative contributions 
of black holes of different final mass, as labelled; the grey area marks the So- uncertainty region of the cosmological star formation rate 
as inferred by Hopkins & Beacom (2006), scaled by a factor of 6.5 X 10~*. Lower panels: predicted cumulative black hole mass density 
as a function of redshift {solid lines), and contributions of black holes of different mass, as labelled. The grey bar indicates the values 
and systematic uncertainties in the total local mass density in black holes estimated by SWM 



lation in bins of mass. In the upper panels, solid lines show 
the growth rate of the integrated black hole mass density 
dpBn/ dt{z) as a function of redshift, while the other lines 
show the mass density accreted in selected bins of current 
black hole mass, as labelled. As already noted by several 
groups (e.g., Marconi et al. 2004; Merloni et al. 2004; Merloni 
& Heinz 2008; SWM), the total pw{z) closely matches the 
shape of the cosmological star formation rate (SFR), here 
taken from Hopkins & Beacom (2006; with the gray area 
marking their S-cr contours) and re-scaled by an ad hoc fac- 
tor of 6.5 X 10"*, close to the ratio between black hole mass 
and stellar mass measured in the local Universe (e.g., Haring 
& Rix 2004). All three models again show a clear signature 
of downsizing in their accretion histories. The accretion onto 
the very massive black holes with \og M^yi/Mq > 8 always 
peaks at z ~ 2, concurrent with the peak in the emissivity 
of luminous optical quasars (e.g., Osmer 1982; Richards et 
al. 2006; Groom et al. 2009). The less massive black holes 
with 7 < log Mbh /Mq < 8 are characterized by a much 
broader peak centered at z ^ 1 — 1.5. The lower panels show 
the cumulative mass density of all black holes with mass 
10^ < Mbh/Mq < 10^° (solid lines), and the mass density 
accreted onto black holes of different current mass, as la- 
belled. By construction, all models share the same radiative 
efficiency and therefore accumulate the same total mass den- 
sities at any time. They differ only slightly with respect to 
the total mass accumulated in different mass bins. The solid 



grey square indicates the systematic uncertainties in the to- 
tal local mass density in black holes estimated by SWM. 

Figures[3]and|4]show that the evolution of the black hole 
mass function is insensitive to the shape of P{X), provided 
the characteristic value Ac and the input AGN luminosity 
function are held fixed. This reassuring result indicates that 
earlier studies assuming a single A value (e.g., Marconi et al. 
2004; Shankar et al. 2004; SWM) reached robust results for 
mass function evolution. This is essentially because in our 
G+P model, with power-law slope a — —0.3, the emissiv- 
ity per logarithmic interval of A is dominated by the Gaus- 
sian peak (Fig. [TJd). A steeper slope would yield duty cycles 
greater than unity and is therefore not allowed. Even our 
G-f P model is disfavored at high redshifts because the ob- 
served quasar clustering implies that the most massive black 
holes should be active at large values of A with high duty 
cycles. Our model results are obviously dependent on the 
value of Ac, which changes the location of the break in the 
black hole mass function, and the value of e, which affects 
the normalization (see SWM and Section |3] below). 



3.3 Comparison with measured P(A) distributions 

Drawing on virial mass estimators grounded in reverberation 
mapping studies (e.g., Wandel, Peterson & Malkan 1999; 
Vestergaard 2002; McLure & Jarvis 2004) and the avail- 
ability of large samples of quasar spectra from wide field 
surveys, several groups have inferred the distribution of Ed- 
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Figure 5. Eddington ratio distributions predicted by the G (solid lines) and G+P (long-dashed lines) models, compared to the K06 
data (grey histograms) for black holes in the range 10^ < AIbh/Mq < lO^*^ at different redshift and luminosity bins, as labelled. The 
predictions have been estimated at the average redshifts of 2 = 1 and 2 = 2 for the low and high-2 subsamples, respectively, although 
the predictions do not strongly depend on these choices if no redshift dependence in the P(A) is included. The cyan long-dashed line in 
the lower-right panel is the G+P model prediction when no lower luminosity threshold is considered. 




Figure 6. Eddington ratio P(X) distributions predicted by the G (solid lines) and G+P (long-dashed lines) models at lower redshifts. 
Left panel: Predictions compared to the H09 data (grey histogram) for black holes in the range 10^ < Mbh/Mq < IQi" with logL > 43.5 
at z = 0.5. Right panel: Predictions compared to the KH09 data (grey histogram) for black holes in the range 10^ < Mbh/A/0 < 
with logL > 42.5 at 2 = 0.2. 



dington ratios in different ranges of redshift and luminosity. 
KoUmeier et al. (2006; hereafter K06) measured P{X) from 
the AGN and Galaxy Evolution Survey (AGES; Kochanek 
et al. 2011), finding that luminous AGNs at 0.5 < 2; < 3.5 
have a quite narrow range of Eddington ratios, with a peak 
at A ~ 0.25 and a dispersion (including observational errors) 
of ~ 0.3 dex. Netzer et al. (2007), analyzing quasar sam- 
ples from the SDSS, found a similar result, with a slightly 
larger dispersion, from a sample centred at 2; '-^ 2.5. Netzer 



& Trakhtenbrot (2007, see also Vestergaard 2004; McLure & 
Dunlop 2004), analyzing nearly ten thousand SDSS quasars 
in the redshift range 0.3 < 2; < 0.75, confirmed the log- 
normal shape of the Eddington ratio distribution and also 
found evidence for a significant decrease with time in the 
characteristic Ac defining the peak of the distribution. K06 
also found a factor ~ 3 decrease in the peak value of A 
between high and low redshifts for the more massive black 
holes in their sample. Many other studies of the Edding- 
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ton ratio distributions of quasars and active galaxies have 
found evidence for time evolution of Ac, and in some cases 
for mass-dependence (e.g., McLure & Dunlop 2004; Vester- 
gaard 2004; Heckman et al. 2005; Ballo et al. 2007; Babic 
et al. 2007; Bundy et al. 2008; Rovilos & Georgantopoulos 
2007; Cao & Li 2008; Fine et al. 2008; Shen et al. 2008; 
Greene et al. 2009; Hickox et al. 2009; Kauffmann & Heck- 
man 2009; Kelly et al. 2010; Shankar et al. 2010e; Stein- 
hardt & Elvis 2010a,b; Trakhtenbrot et al. 2011; Willott et 
al. 2010a,b; Aird et al. 2011; Shen & Kelly 2011). However, 
sample selection, observational noise, and intrinsic scatter in 
black hole mass estimators can all have strong effects on the 
apparent distribution of Eddington ratios (e.g., Lamastra et 
al. 2006; Lauer et al. 2007; Marconi et al. 2008; Shen et al. 
2008; Netzer 2009; Shen & Kelly 2010; Rafiee & HaU 2011). 
At low redshifts, where samples can probe to much lower 
AGN luminosities, it is clear that the distribution of Ed- 
dington ratios is much broader than the roughly log-normal 
distribution found for optically luminous quasars at high 
redshift. 

As a set of representative examples of the literature, we 
consider the KoUmeier et al. (2006; K06 hereafter) results 
for broad-line optically luminous quasars aX z > 0.5, the 
Hickox et al. (2009; H09 hereafter) results at z ~ 0.5 for a 
sample selected by X-ray, infrared, and radio emission, and 
Kauffmann & Heckman's (2009; KH09 hereafter) analysis 
of the local AGN population in the SDSS, with an effective 
redshift Z7»G.2. 

To make a close comparison of our models to observed 
samples of AGNs, we compute Eddington ratio distributions 
conditioned to active AGNs with luminosities between Lmin 
and i^max , -Pobs , as 

-Pobs(A|2, -£/niin, -£/max) = (13) 

^ X] X] P(A = Lk/{MmMz, MBH.j)iVact (Mbhj, z) , 

k j 

where N^ctiMB-a, z) is defined in Eq. [21 The sums in the 
above equation are extended to all bins of black hole mass, 
and to all bins of luminosity in the assumed observed range, 
L-mi-n < Lk < I/max- Siucc we are interested here in the 
shapes of the predicted Eddington distributions, the con- 
stants A are chosen to renormalize the Pobs(A|2:, I/min, imax) 
to match the peaks of the observed distributions. 

Fig. [S] compares the Eddington ratio distributions pre- 
dicted by the G and G-l-P models to the K06 data, for black 
holes in the range lO'^ < Mbh/Mq < 10^° in four different 
bins of redshift and luminosity. The predictions have been 
computed at the redshifts of 2; = 1 and 2 = 2 for compari- 
son to the z < 1.2 and z > 1.2 subsamples, respectively. The 
predictions do not strongly depend on these choices, since 
P(A) has no explicit redshift dependence. Both models pre- 
dict, essentially by construction, a roughly log-normal P(A) 
with a typical Ac that agrees with the K06 measurements. 
Most importantly, the power-law tail of the G-l-P model does 
not lead to substantial disagreement because the fraction of 
quasars of a fixed luminosity with low values of A is much less 
than the fraction of black holes of a fixed mass with the same 
low values of A. In fact, the A distribution for all quasars with 
any luminosity (cyan dashed curve in the lower right panel) 
has the same shape as our assumed intrinsic distribution in 
Figure la. The most significant disagreement with the data 



is the high characteristic A for 45.5 < logL < 46 at 2 < 1.2 
(lower left panel), a consequence of keeping the model Ac 
fixed instead of allowing redshift evolution. 

Fig. |S] compares the predicted and observed A distribu- 
tions for the lower redshift samples of H09 at 2 ~ 0.5 and 
KH09 at 2 ~ 0.2. To roughly take into account the efi'ective 
limits in black hole mass and luminosity of the observational 
samples, for the former we consider all black holes with 
Mbh > lO'' Mq shining at L > 10'*^ergs"\ while in the lat- 
ter case we adopt the restricted range 10^ < Mbu/Mq < 10® 
shining above L — lO''^'^ ergs~^. The match to the observa- 
tions is poor in both cases. The distributions predicted by 
the G model are much narrower than the observed distribu- 
tions, and they peak at a value Ac ~ 0.25 that is roughly ten 
times higher than the high-A peaks of the observed distri- 
butions. The power-law tail in the G-l-P model is in rough 
agreement with the observed distributions at low A, but the 
peak of the Gaussian component remains discrepant. Fitting 
these observations, even approximately, requires evolution of 
Ac, decreasing towards low redshift. We return to this point 
in Section [5. 21 b elow . after first finding independent evidence 
for evolution in Ac from observed active galaxy fractions. 



4 ACTIVE GALAXY FRACTIONS: REDSHIFT 
AND MASS-DEPENDENCE OF P(A) 

Duty cycles are a basic prediction of continuity equation 
models as discussed in Section [2l If we assume that most or 
all massive galaxies contain black holes, as supported by ob- 
servations of the local universe (see, e.g., Ferrarese & Ford 
2005, and references therein), then active galaxy fractions 
can be used as an observable proxy for black hole duty cy- 
cles. However, active galaxy fractions depend on the lumi- 
nosity or Eddington ratio threshold of each specific observa- 
tional sample: lower thresholds obviously increase the active 
fraction. Therefore, in order to properly compare model pre- 
dictions to available data sets, we define the fraction of black 
holes shining above a luminosity threshold Lmin , or observed 
duty cycle l7obs, as 

/"OO 

t/ob.(AfBH,Lmin,2) = / P {X\MbH , z)U {MbH , z)dlog X , (14) 

J log Ai-i-^in 

Amin ~ Lmin/ilMBu) ■ 

In this Section and Section [5] we will introduce versions 
of our reference models that incorporate redshift and mass 
dependence of P(A), and we will consider a variety of "test" 
models that illustrate specific points. The key features of 
these models are summarized in Table [TJ and we will sum- 
marize the qualitative successes and failures of the reference 
models in Table [2] below (Section ??). These model changes 
generally have little impact on the overall mass accretion his- 
tories, so we will not repeat the analysis shown previously 
in Figures [3] and [4] but instead focus on predictions where 
the new models differ significantly from those described in 
Section [31 Despite systematic uncertainties in the observa- 
tional constraints, even qualitative trends with redshift and 
black hole mass are enough to provide strong model tests. 
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Figure 7. Comparison of predicted and observationally estimated duty cycles at low redshift. Left panel: Data points from Goulding 
et al. (2010) compared to model predictions at z = 0.2 as labelled, with a luminosity tfireshold logL > 41.5. Middle panel: Data points 
from Kauffmann et al. (2003) compared to model predictions, with a luminosity threshold logL > 43.5. Right panel: Data points from 
Schulze & Wisotzki (2009) compared to model predictions, with an Eddington ratio (not luminosity) threshold A > 0.01 and a correction 
factor for the Type 1 Broad Line AGNs taken from Greene & Ho (2009). The three lowest mass points include the Schulze &; Wisotzki 
(2010) incompleteness corrections, which are negligible at higher mass. 



4.1 Observational Estimates 

We first summarize the observational estimates of AGN frac- 
tions that we adopt for our model comparisons, with fur- 
ther details of some of these estimates given in Appendix [D] 
For high luminosity, broad-line (Type 1) AGN in the lo- 
cal Universe, Schulze & Wisotzki (2010) have estimated the 
mass function of active black holes at z < 0.3 by applying 
linewidth mass estimators to quasars identified in the Ham- 
burg/ESO objective prism survey. They divide their active 
black hole density by the total black hole space density in 
the Marconi et al. (2004) mass function to derive an active 
fraction. Schulze & Wisotzki impose a threshold of A ^ 0.01 
in their estimated Eddington ratio to define active systems. 
Host galaxy contamination causes incompleteness in their 
AGN catalog below A/bh ~ lO'' '' Mq (for A ^ 0.01). In the 
right panel of Fig. [7] we plot the points from their Fig. 12, 
correcting their three lowest Mbh points for incompleteness 
as suggested by the authors. According to their simulations, 
higher Mbh points should be unaffected by incompleteness, 
and we have omitted points below Mbh = lO^M© where 
incompleteness corrections become large. The active black 
hole fractions estimated by Schulze & Wisotzki (2010) are 
low, only ~ 2 X 10"" at A/bh = l(f Mq, rising to ~ 10"^ at 
Mbh = 10* Mq. 

Kauffmann et al. (2003) estimate much higher active 
fractions for narrow line (Type 2) AGN identified by apply- 
ing emission line diagnostics to SDSS spectra in the SDSS 
main galaxy redshift survey. Since Kauffmann et al. (2003) 
begin with a galaxy catalog, quasars that are bright enough 
to appear as point sources (rather than extended sources) 
in SDSS imaging are excluded. They estimate their com- 
pleteness limit at 10^ Lq in [OIII] luminosity, which we con- 
vert to an approximate Lmin ~ lO^^'^ergs"^ by adopting an 
extinction-corrected bolometric correction of ~ 800 (both 
Kauffmann et al. 2003 and KH09 adopt extinction-corrected 
luminosities). Points in the middle panel of Fig. [7] are taken 
from Fig. 5 of Kauffmann et al. (2003), where we have sim- 
ply converted stellar mass to black hole mass by multiplying 
by 1.6 X 10~^ (Magorrian et al. 1998). This crude conversion 
should be reasonably accurate above Af, ~ 10^^'^ A/0, but it 
probably overestimates the black hole mass for lower mass 



galaxies that are no longer bulge-dominated. The Kauff- 
mann et al. (2003) active fractions rise from ~ 2% at 
A/bh = 10^ M© to ~ 10% at Mbh = lO^Mg, then decline 
towards lower A/bh. The 10'*"^'^ergs~^ luminosity threshold 
corresponds to A = 0.01 at A/bh ~ 10^'^ M©, so at higher 
masses the Kauffmann et al. (2003) estimates correspond to 
a lower A-threshold than Schulze & Wisotzki's. 

The 1.5 to 2 order-of-magnitude gap between the Kauff- 
mann et al. (2003) and Schulze & Wisotzki (2009) active 
fractions can be partly explained by the different abun- 
dances of broad line and narrow line AGN at low redshift. 
The abundance ratio of the two types of AGN has been ad- 
dressed by Greene & Ho (2009; a correction of Greene & 
Ho 2007), who analyze all SDSS spectroscopic objects (tar- 
geted as galaxies or quasars) in Data Release 4 (Adelman- 
McCarthy et al. 2006) to define a sample of 8,400 broad-line 
AGN at 2; < 0.3. Their estimated active fractions (in their 
Fig. 11) are close to those of Schulze & Wisotzki. They find a 
gap of 0.9-1.3 dex between their active black hole mass func- 
tion for broad-line AGN and the corresponding result (Heck- 
man et al. 2004) for SDSS narrow-line AGN (see Greene & 
Ho 2009, Figure 10), which implies a correction of ~ 10 — 20 
for the ratio of Type I to Type II AGN at low redshift. 
The remaining factor of 5 — 10 may be accounted for by the 
effects of scatter in bolometric corrections, different luminos- 
ity thresholds, and systematic offsets in the black hole mass 
function estimates. Note that scatter in /'boi//'[oni] could 
make the effective luminosity threshold of the Kauffmann et 
al. (2003) data lower than we have assumed, in particular if 
the luminosity function is affected by a steeply rising low-A 
tail in P(A). 

For studies that probe to lower AGN luminosities, active 
galaxy fractions are even higher than those of Kauffmann et 
al. (2003). For example. Ho (2004) argues that at least 40% 
of local bulge-dominated galaxies host a low luminosity AGN 
and/or a LINER, and Grier et al. (2011) find nuclear X-ray 
sources in ~ 60% of galaxies in the SINGS survey (Kennicutt 
et al. 2003). As a representative example of low luminosity 
statistical studies, we show in the left panel of Fig. [7] the 
data of Goulding et al. (2010, their Fig. 5), based on an 
X-ray and IR census of a volume-limited sample of galaxies 
with D < 15 Mpc. Their active fractions range from 14% to 
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Figure 8. Predicted and observationally estimated duty cycles at redshifts z = 0.5, z = 1.3, and z = 2.3. In all panels, grey bands show 
the estimates of Xue et al. (2010). Points in the middle panel, points are from Bundy et al. (2008). Points in the right panel come from 
Erb et al. (2006; filled small squares), Kriek et al. (2007; filled triangles), Alexander et al. (2005; 2008; filled circle), and Caputi et al. 
(2006; semi-filled circle). 



37% over the black hole mass range 10^ - lO*'* Mq, albeit 
with large statistical errors reflecting the small sample size. 
The luminosity threshold of this sample is very low, about 
LuAtl = 10*^'^ ergs~^ in bolometric units (their Fig. 3). 

At higher redshifts we take as our primary data set the 
measurements reported by Xue et al. (2010) based on mod- 
erate luminosity X-ray AGN in the Chandra Deep Fields. 
We reproduce results from their Fig. 14 (the upper, orange 
bands) as grey bands in the three panels of Fig. (8] convert- 
ing their X-ray luminosity thresholds to bolometric lumi- 
nosity thresholds with our luminosity-dependent bolomet- 
ric correction (Section l2.2|) . We again convert stellar masses 
to black hole masses by simply multiplying by the Magor- 
rian et al. (1998) factor of 1.6 x 10^^. This scaling should 
be taken with a grain of salt, as we are using total stel- 
lar masses rather than bulge masses and ignoring possible 
redshift evolution of the scaling factor. However, these un- 
certainties should not affect the key lessons that we take 
from Xue et al. (2010): at all three redshifts, the duty cycle 
for BHs shining above logL > 43.3 increases with increasing 
Mbh (indicated in the data by increasing M«), and for the 
most massive galaxies active fractions are 10 — 50% (grey 
bands marking the 1-a uncertainty). 

Other points in Fig. |8] come from other studies with 
similar redshift ranges and luminosity thresholds. At in- 
termediate redshifts, the data of Bundy et al. (2008) from 
the DEEP2 and AEGIS surveys show nearly the same mass 
trend and normalization as Xue et al. (2010). More recently, 
Mainieri et al. (2011) also found evidence for a steep in- 
crease in AGN activity with stellar mass at (z) ~ 1.1, 
though with a luminosity threshold about two orders of 
magnitude higher than the one by Xue et al. (2010). At 
« > 2 we add several points for high mass black holes based 
on the studies of Caputi et al. (2006), Erb et al. (2006), 
Kriek et al. (2007), and Alexander et al. (2005, 2008). These 
data have been collected at bolometric luminosities approx- 
imately log L > 43.8, comparable to those of Xue et al. 
(2010). Further details about these observations axe dis- 
cussed in Appendix|D| They extend the mass trend found by 
Xue et al. (2010), with active fractions ranging from ~ 10% 
to ~ 50%. 



4.2 Low-redshift AGN fractions: A 
redshift-dependent P(A)? 

Red dotted and cyan dot-dashed curves in Fig. [7] show the 
duty cycles C/(Mbh, z) predicted by the G and G-l-P models 
at 2 = 0.2, where we have imposed the bolometric luminosity 
thresholds logL > 41.5 (left panel), logL > 43.5 (middle), 
and A > 0.01 (right). For the right panel, we have addi- 
tionally multiplied model predictions by a mass-dependent 
factor (ranging from ~ 15% at Mbh = IO^Mq to ~ 3% at 
Mbh = 3 X 10** M©) to account for Greene & Ho's (2009) 
estimate of broad line AGN fractions. (Specifically, the cor- 
rection is obtained from the ratio of Broad Line to Type 
II number densities in their Figure 10, with an extrapola- 
tion of their results above A^bh ~ 5 x IO^Mq.) Both mod- 
els drastically underpredict all measured active fractions for 
Mbh > 10^' ^Mq. Predicted duty cycles are higher for the 
G-l-P model, and the gap between G-l-P and G is larger for 
lower luminosity thresholds that allow a larger contribution 
of sub-Eddington black holes, but the G-l-P model still falls 
well short of the observed active fractions. 

As introduced by Shankar et al. (2004) and discussed 
more extensively by SWM (see their Figure 11), one can in- 
crease predicted duty cycles at low redshift by adopting a 
redshift-dependent Eddington ratio distribution, with char- 
acteristic Ac that drops towards low redshifts. Decreasing 
Ac oc L/Mbu maps a given luminosity to more massive, and 
hence rarer, black holes, thus requiring a higher duty cy- 
cle to reproduce the observed AGN space density. Further- 
more, the results of H09, KH09, and other studies directly 
suggest a lower Ac at low redshifts, as already discussed in 
Section r3.3l (see Fig.|6]). Motivated by these results, we intro- 
duce models in which the location of the peak in the G and 
G-fP -P(A) distributions decreases from Ac = 1.0 at z = 6 
to Ac ~ 0.02 at 2 0.1, following 



Ac(2) = A(2 = 6) 



1 + 2 



7 



(15) 



with a = 2.2, although we note that even lower values of 
a yield similar results. Eg. ITSl implies Ac = 0.71, 0.48, 0.29, 
0.16, 0.06, 0.03, 0.02 at 2 = 5, 4, 3, 2, 1, 0.5, and 0.1, re- 
spectively. 

In agreement with SWM, who considered (5-function 
P(A), we find that this evolution has a minor impact on the 
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evolved black hole mass function (at least at high masses) or 
on the global accretion histories shown in Fig. U However, 
predicted duty cycles at 2 ^ 0.2 for the G(z) model — a 
Gaussian P(A) with the Xc{z) given by Eq. [15] — are much 
higher than those of the G model, as shown by comparing 
the dotted and long-dashed curves in Fig. [71 The model now 
approximately agrees with the Goulding et al. (2010), Kauff- 
mann et al. (2003), and Schulze & Wisotzki (2010) active 
fractions, though at high masses it is above the latter data 
and below the former two. Solid and short-dashed curves in 
Fig. [7] show models with Gaussian and G-l-P A-distributions 
that incorporate both redshift and mass dependence of Ac, 
as discussed in the next section. 

The redshift dependence expressed in Eq. [15] is by no 
means unique. In particular, we have checked that an equally 
good match to the local data can be obtained by allowing the 
characteristic Ac to decrease only below redshift 2 ~ 1, while 
remaining close to unity at higher redshifts. This abrupt 
change in Ac at late time^ leads to duty cycles that increase 
at fixed black hole mass from 2 ~ 1 to 2 ~ 0, instead of de- 
creasing as predicted by continuous redshift evolution. This 
increasing behaviour may be in conflict with observations; 
some direct X-ray and optical analyses show that the frac- 
tion of active galaxies above a fixed stellar mass decreases 
with decreasing redshift in the range < 2 < 1 (e.g., Shi et 
al. 2008). 

4.3 AGN fractions vs. mass: A mass-dependent 

P(A)? 

Fig. [8] compares the model predictions to the higher redshift 
active fraction data discussed in Section 14.11 The G and 
G(2) models both predict a duty cycle that declines with 
increasing A/bh, opposite to the trend found by Xue et al. 
(2010) and Bundy et al. (2008). At first glance, the G+P 
model appears to fare better, at least in the 2 = 1.3 and 
2 = 2.3 panels. However, we view this better agreement as 
artificial — it is a consequence of the high initial black hole 
space density required to keep high-2 duty cycles below unity 
in this model (see Fig.[3|. The high black hole space density 
inherited from these initial conditions leads to a lower duty 
cycle for low mass black holes. We regard these high initial 
space densities as observationally untenable for the reasons 
discussed in Section[3]2] so we discount this apparent success 
of the G+P model, which vanishes by 2 = 0.5 in any case. 

The only other way we have found to make the predicted 
duty cycle increase with black hole mass at these redshifts 
is to introduce an explicit mass dependence of -P(A), with 
Ac declining with increasing mass. In this case massive black 

^ We note that Shankar, Bernardi & Haiman (2009) found that 
a constant Eddington ratio as a function of redshift was consis- 
tent with the statistics of local, early-type galaxies coupled with 
a mild evolution in the black hole-velocity dispersion scaling re- 
lation and the integrated AGN energy density. However, their 
conclusions were driven by assuming a non-evolving structural 
evolution of their hosts. Allowing for some evolution in velocity 
dispersion (as suggested by several observations and models, as 
in, e.g., Shankar et al. 2011, and references therein) would be 
consistent with a decrease in time of the characteristic Eddington 
ratio while preserving the mild evolution in the black hole- velocity 
dispersion relation. 



holes are matched to more common, lower luminosity AGN, 
implying a higher duty cycle. There is some direct empirical 
support for such a trend (e.g. Heckman et al. 2004; Ballo 
et al. 2007; Babic et al. 2007; Netzer & Trakhtenbrot 2007; 
Rovilos & Georgantopoulos 2007; Fine et al. 2008; H09), and 
it could arise theoretically in models that envisage a faster 
shut-off of activity in more massive systems due to AGN 
feedback (e.g., Matteucci 1994; Granato et al. 2006; Lapi et 
al. 2006). Motivated roughly by these empirical studies, we 
have adopted a model with the same redshift dependence as 
before and a mass dependence 

log Ac(Mbh, 2) = log Ac(2) + [/3M(log Mbh - C)] , (16) 

with Pm — —0.3 and C = 7. Eg. 1 16 1 implies a Ac that drops 
by a factor of two between IO'^Mq and 10* M© and by a 
further factor of two between 10* Mq and 10^ Mq. 

Solid lines in Figures [7] and [8] show predictions of the 
G(2,A'/bh) model, which incorporates mass-dependence of 
the radiative efficiency e as well as of Ac, for reasons that we 
will discuss shortly. This model also incorporates a broad- 
ening of the Gaussian in the G(2,Mbh) at lower redshifts 
following the trend 

logEA(2) = logEA(2 = 6)-0.41og(^i±^^ , (17) 

with logEA(2 = 6) = 0.3, yielding T.x(z) ~ 0.4,0.5,0.6 at 
2 = 3, 2, 0.2. As discussed below, we include this addi- 
tional modification to the model to provide a better match 
to the local Eddington ratio distributions. At all redshifts, 
this model predicts a duty cycle that is roughly fiat or ris- 
ing with black hole mass, though the rising trend is always 
weaker than that found by Xue et al. (2010) and Bundy et 
al. (2008). We could adopt a still stronger mass dependence 
in this model, which would improve agreement with the ob- 
served trends, but this would no longer be supported by 
empirical estimates (e.g., Hickox et al. 2009), and it would 
exacerbate problems in explaining the Eddington ratio dis- 
tributions and local black hole mass function (see below). 
The G(2,AfBH) model is also in rough agreement with all 
the local data sets on active fractions (solid lines in Fig. [7]), 
though it overpredicts the Schulze & Wisotzki (2010) active 
fractions for Mbh ~ 10^ Mq by a factor of several. 

Short-dashed lines in Fig. [H] show the predictions of a 
model labeled G-fP(2,MBH) that has a G+P form with 
redshift- and mass-dependent Ac and mass-dependent e like 
that of the G(2, A/bh) model. Instead of allowing the Gaus- 
sian to broaden at low redshift via Eq. 1171 we keep its width 
fixed at Ea = 0.3 but allow the power-law tail to grow over 
time. Specifically, the P(A) distribution is cut off at 

A„,in(2) = 0.1 (^i±^y , (18) 

implying A^in = 0.1, 0.063, 0.036, 0.019, 0.0079, 0.0023, 
0.0098, 0.00039 at 2 = 6, 5, 4, 3, 2, 1, 0.5, 0.1. This simple 
modification (similar to that adopted by Cao 2010) avoids 
the need for a large initial black hole seed population, thus 
removing the artificial aspect of our redshift-independent 
G+P model, but it produces a full G+P distribution at in- 
termediate and low redshifts. The predictions of this model 
are similar to those of G(2,AfBH), with a somewhat worse 
match to the 2 = 1.3 data in Fig. [Sjj and a somewhat better 
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Figure 9. P(A) distributions for our G{2, A/bh) model {dashed curves) and the G(2,A/bh) model with a steep power-law component 
added at low redshift (solid curves), shown at z = 0.6 {left) and z = 0.2 (right). In each panel, the three curves show P(A) for black hole 
masses of IO'^Mq, IO^Mq, and IO^Mq from top to bottom, as labelled. 



match to the Schulze & Wisotzki data in Fig. [T] Both mod- 
els fail to reproduce the rising U{Mbh) found at z = 0.5 by 
Xue et al. (2010). 

4.4 AGN fractions vs. mass: a steep power-law 
component at low z? 

At early stages of our investigation, we anticipated that the 
observed trends of rising active fraction with rising galaxy 
mass might be explained mainly by the combination of a 
broad P(A) with sample luminosity thresholds: high mass 
black holes remain observable at lower A values, so they 
could have higher duty cycles above fixed luminosity even if 
the trend of total duty cycle (over all A values) was flat or 
decreasing. This is essentially the explanation advanced by 
Aird et al. (2011), who fit AGN data from the PRIMUS sur- 
vey with a model in which all black holes have an Eddington 
ratio distribution P(A) = Po{z) x (A/Ao)~°'^ independent of 
mass. 

As already discussed in the context of our G+P model, 
we find that a long power-law tail is unrealistic at high red- 
shift because it implies duty cycles higher than unity unless 
the population of seed black holes is implausibly large. How- 
ever, a power-law component that kicks in at lower redshifts 
is still feasible. The blue dotted lines in Figures [7] and [8] 
show the predicted duty cycles for a model equivalent to 
G{z, Mbh) with the addition of a power-law component at 
z < 0.7 with slope a = —0.9 and normalized to have the 
same value as the Gaussian at log A = log Ac — 0.2 log E^. 
Fig. [5] plots P(A) for this model at three different black 
hole masses and two redshifts, in comparison to those of 
the G{z, Mbh) model. For this plot we have multiplied P(A) 
by U{Mbh, z) so that the curves are normalized to the duty 
cycle rather than integrating to unit probability. The Gaus- 
sian and power-law components join to produce a P{\) that 
is close to a single power law from A — to A = 1. 

This model still exhibits "downsizing" in the sense that 
higher mass black holes have lower duty cycle at any given 
A. However, the duty cycle for AGN active above a lumi- 
nosity threshold logL = 43 increases with black hole mass, 
as shown in the left panel of Fig. |8l because massive black 



holes can shine above the threshold at low A. This is the only 
model we have constructed that reproduces the Xue et al. 
(2010) trend at 2 = 0.5. The Aird et al. (2011) prescription, 
with the P{\) normalization independent of mass, would 
produce a still stronger rising trend, but we find that such a 
model cannot simultaneously match our evolved black hole 
mass function and our input luminosity function. 

5 OTHER OBSERVATIONAL CONSTRAINTS 

5.1 The local mass function: A mass-dependent 
radiative efficiency? 

As shown previously in Fig. O our simple, non-evolving 
5, G, and G-l-P models yield reasonable agreement with 
SWM's estimate of the local black hole mass function for 
Mbh ^ 10^ Mq, but the model predictions skirt the upper 
boundary of the observational estimates at higher masses. 
Fig. [To] plots mass functions for the Gaussian P(A) mod- 
els with various assumptions about redshift evolution and 
mass-dependence. Redshift evolution alone (red long-dashed 
line) produces modest changes above 2 x 10* Mq, but the 
decreasing Ac boosts growth of high mass black holes at 
the expense of low mass black holes, leading to a flatten- 
ing of $(Mbh) near 10** and a lower space density of 
lower mass black holes. This change partly reproduces the 
turnover found by Vika et al. (2009) , though it still does not 
reproduce their estimate below Mbh ~ lO^'^M©. 

Adding only the mass- dependent Ac of Eq. [16] leads to 
the "G(2:, Mbh), constant e" curve (blue dot-dashed line) in 
Fig. 1101 For this model, we have lowered the radiative ef- 
ficiency from e = 0.06 to e = 0.05 to improve the match 
to the amplitude of $(Mbh) at 10* - IQ^'^Mq, where it is 
best measured. However, the model then strongly overpre- 
dicts the mass function at A/bh > 10^ Af©, a direct conse- 
quence of the higher duty cycles of high mass black holes 
that the model was designed to produce. As discussed in 
Section 15.31 below, the expected impact of black hole merg- 
ers would make this overprediction even more severe. 

Volonteri, Sikora & Lasota (2007), Cao & Li (2008), and 
Fanidakis et al. (2011), among others, have suggested that 
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Figure 10. Local black hole mass function predicted by different models, as labelled. Data are as in Fig. [3] Of the three G(2,Mbh) 
models, with mass-dependent Ac, only the one with mass-dependent e (solid line) matches the high end of the mass function. Elsewhere 
in the paper this preferred model is simply labeled G{z, Mbh)- 



radiative efficiency may increase with black hole mass, per- 
haps mirroring a merger-induced increase in the average spin 
for the more massive black holes. Cao & Li (2008) claimed 
evidence for an increasing e from the match between the pre- 
dicted and local black hole mass functions at the high-mass 
end. More recently Davis & Laor (2011) directly determined 
the radiative efficiency from the ratio between bolometric lu- 
minosity and accretion rate, the latter determined from thin 
accretion disk model fits to the optical luminosity density. 
Their analysis seems to support an increase of the radiative 
efficiency with black hole mass, from 0.03 at low masses to 
0.4 at high masses. An increasing radiative efficiency at high 
masses is also consistent with the notion that the bulk of lu- 
minous radio AGNs, believed to be rapidly spinning black 
holes, are indeed massive black holes (e.g., McLure & Jarvis 
2004; Metcalf & Magliocchetti 2006; Shankar et al. 2008a,b; 
Shankar et al. 2010e). 

Inspired by these theoretical and empirical arguments, 
we have constructed the G(2:,Mbh) model (see Tables [1] 
and [21, which, following Cao & Li (2008), adopts a mass- 
dependent e given by: 

^ r 0.05 if Mbh < W^Mq 

" \ O.O5(Mbh/1O**M0)°-3 ifMBH^lO^Me . ^ ^ 

This change reduces the implied accretion rates of massive 
black holes, leading to much better agreement with the local 
$(Mbh) as shown by the solid curve in Fig. 1101 This model 
also has a stronger positive trend of duty cycle with Mbh, 
improving agreement with the data in Fig. [8] The two im- 



provements are connected: raising e decreases the space den- 
sity of massive black holes, so a higher duty cycle is required 
to match the space density of active black holes. The mass 
function (not shown) of the model introduced in Section [4.41 
with a steep power-law at 2 < 0.7, is nearly identical to that 
of G(z,Mbh), except at log A/ ^ 7.2, where it is lower by 
~ 0.3 dex. 

Wang et al. (2009) have instead recently claimed em- 
pirical evidence for an increase of e with redshift. Their 
results are based on an inversion of Soltan's (1982) argu- 
ment, expressing the radiative efficiency at any redshift as 
the ratio of the accreted mass density up to that redshift 
Pbh(z) to the corresponding total emissivity obtained by 
direct integration of the AGN luminosity function (see their 
Eq. 6). The mass density pbh(2;) was computed by first mea- 
suring the mass density locked up in all active black holes 
at z (with masses from virial relations), then correcting by 
the duty cycle extracted from the number counts of active 
galaxies in VIMOS-VLT Deep Surveys. Clustering analysis 
also provides hints of redshift-dependent radiative efficiency. 
Shankar et al. (2010b), adopting basic accretion models and 
cumulative number matching arguments, found that black 
hole accretion plus merger models consistent with both the 
quasar luminosity function and the strong observed cluster- 
ing at 2 « 4 (Shen et al. 2007) must be characterized by high 
duty cycles and large radiative efficiencies e > 0.2, if they 
are accreting at a significant fraction of the Eddington limit. 
However, an efficiency e ^ 0.2 at all redshifts would under- 
predict the local black hole mass function (SWM). In fiux- 
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limited quasar surveys, higher redshift quasars are found 
above higher luminosity thresholds, so both of these em- 
pirical arguments could potentially be answered by a mass- 
dependence of the sort implied by Eq. 1161 

We nonetheless consider a directly redshift-dependent 
model, approximating these empirical findings with the re- 
lation 



e{z) = 0.0022 



1 -l-erfc - — 



z 



(20) 



Eq. (O implies e ^ 0.14, 0.14, 0.12, 0.05, 0.02, 0.01 at 
2 = 6, 4, 2, 1, 0.5, 0.2, respectively. The local black hole 
mass function predicted by this G{z,Mbw) + model, 
shown by a green dashed curve in Fig. IIOI tends to overes- 
timate the observed mass function at all scales. The match 
could be improved by increasing the overall normalization of 
e[z), still allowed (and actually preferred) by the measure- 
ments of Wang et al. (2009; see their Figure 2a), but such an 
increase produces unphysical models with duty cycles signif- 
icantly higher than unity, and it still leaves too many high 
mass black holes. Moreover, the cosmological accretion rate 
predicted by the G(2:, Mbh) + ^{z) model is morphologically 
different from the cosmological star formation rate of galax- 
ies, at variance with the agreement shown in Fig. 2] In itself 
this is not a fatal objection, but it then requires non-trivial 
fine tuning to reproduce a tight local relation between black 
hole mass and stellar mass. 



5.2 Eddington Ratio Distributions, Revisited 

In any model with a broad -P(A), the observable distribu- 
tion of Eddington ratios depends on the luminosity of the 
AGN being considered. Fig. [11] shows the contribution to 
the bolometric luminosity function from different ranges of 
\ At z = 3.2 (top panels) and z = 0.3 (bottom panels) in the 
G(z), G(z,Mbh), and G+P(^,Mbh) models (left, middle, 
and right panels, respectively). There is a natural trend for 
the high end of the luminosity function (log L > 46) to be 
mainly contributed by high-A black holes, with a broader 
range of A contributing at lower luminosities. Interestingly, 
we find that the A < 0.01 range does not dominate at any 
luminosity at any redshift. Only in the G(2:,AfBH) model 
and at low redshifts are there roughly equal contributions 
to the luminosity function from each logarithmic bin of A 
in the range 0.01 — 1.0. This is important to compare with 
models that rely on a significant contribution from "ADAF- 
type" modes to the global accretion history of black holes 
(see, e.g., Merloni & Heinz 2008 and Draper & Ballantyne 
2010). 

Fig. [T2] compares the Eddington ratio distributions 
of the G(z,A/bh), G(2:), and G+P(2:,Mbh) models (black 
solid, red long-dashed, and cyan dashed lines, respectively) 
to the observational estimates of K06, H09, and KH09, 
discussed earlier in Section 13.31 Beginning at low redshift 
(bottom panels), we see that the declining \c{z) in the 
redshift-dependent models resolves the discrepancy seen in 
Fig. [6l producing much better agreement with the H09 
and KH09 data. The additional broadening in the Gaus- 
sian component (Eq. I17[) produces excellent agreement with 
the II09 distribution and excellent agreement with KH09 
for A > 10~^'^. The G-|-P(2, A^bh) model achieves similar 



agreement, matching KH09 slightly better at low luminosi- 
ties but still falling short at A < 10~^. The blue dotted 
curves in the lower panels show the G(2,A/bh) model with 
the steep power-law at low z, introduced in Section [4.41 The 
predicted Eddington ratio agrees poorly with both the H09 
and KH09 data. 

At 2 > 1.2 (upper panels) model predictions are in ap- 
proximate agreement the K06 data, with G(z) showing the 
best agreement. However, the z = \ model outputs tend to 
disagree with the K06, z < 1.2 histograms (middle panels). 
We note that Netzer & Trakhtenbrot (2007) find a peak at 
A « 0.1 for M = 10^- 10^-^ Mq black holes at z = 0.7. Kelly 
et al. (2010) have also recently claimed an Eddington ratio 
distribution from SDSS of Broad Line Quasars that peaks 
at I//Z/Edd '^-^ 0.05 with a dispersion of ^ 0.4 dex. Both these 
results would be closer to the model predictions. 

While the duty cycle data considered in Section 14.31 
seem to favor a Ac that decreases at higher Mbh, the 
K06 data disfavor this solution. We note that the lumi- 
nosity threshold of the Xue et al. (2010) study that mo- 
tivates this model is much lower than that of the K06 data, 
Ln,in ~ 10*^ erg s"^ rather than Lmin ~ lO"^ - lO'^'erg 8"^ A 
reconciliation of these results could therefore lie in a model 
that behaves differently in these two luminosity regimes. 



5.3 The Impact Of Mergers 

It is usually assumed that mergers of galaxies are followed by 
mergers of their central black holes, making black hole merg- 
ers a potentially important source of evolution in the black 
hole mass function (e.g., Hughes & Blandford 2003; Wyithe 
& Loeb 2003; Islam et al. 2004; Scannapieco & Oh 2004; 
Yoo & Miralda-Escude 2004; Volonteri et al. 2005; Lapi et 
al. 2006; Yoo et al. 2007; MaruUi et al. 2008; Bonoh et al. 
2009; Shen 2009; SWM; Bonoh et al. 2010; Shankar 2010; 
Shankar et al. 2010a,d; Shankar et al. 2011; Kocsis & Sesana 
2011; Kulkarni & Loeb 2011). We will save a full discussion 
of mergers for future work, but here we briefly assess their 
potential impact on the mass function. Improving on our 
simplified calculation in SWM, we here follow the schemes 
proposed by Shen (2009) and Shankar et al. (2010b), com- 
puting the rate of black hole mergers from the halo merger 
rate predicted from fits to N-body dark matter simulations 
(Fakhouri & Ma 2008), corrected by a dynamical friction 
timescale. In the presence of mergers the continuity equa- 
tion reads as 

dUBH,.. ,-. a((AfBH)nBH(A/BH,t)) , „ „ 

- (Mbh , I) = ^777 h Oin — Oou 



dt 



(21) 

where Sin and Sout are, respectively, the merger rate of 
smaller mass black holes ending up with mass A/bh and 
the merger rate of black holes with initial mass A/bh merg- 
ing into more massive systems. The merger rate of haloes 
as a function of redshift, mass ratio of the progenitors, and 
mass of the remnants are taken from Fakhouri & Ma (2008) . 
We then convert the merger rate of haloes to a merger rate 
of black holes via the median A/aH-Aifhaio relation, defined, 
at all times, by cumulative number matching between the 
black hole and halo mass functions, which allows us to as- 
sociate the proper halo merger rate to a given bin of black 
hole mass. Full details are given in Appendices |B] and [C] 
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Figure 11. Contributions of different Eddington ratio distributions to the overall bolometric luminosity functions at z = 3.2 (top) and 
2 = 0.3 (bottom), as predicted by the G(2), G(2,Mbh)> S'lid G+P(2,Mbh) models (left, middle, right panels, respectively). 



We make the limiting case assumption that halo mergers 
are always followed by black hole mergers after a dynami- 
cal friction time. The accuracy of this assumption remains 
a matter of debate (e.g., Cavaliere & Vittorini 2000; Shen 
2009; Shankar 2010, and references therein). The impact of 
the dynamical friction time delay itself is irrelevant at z <2 
(see Shen 2009). 

Fig. ll3l shows the mass function of the G(z, Mbh) model 
at z — 4, 2, and 0.1 without mergers (dashed lines), includ- 
ing major mergers above a black hole mass-ratio threshold 
5 — 0.5 (solid lines), and including all mergers above a black 
hole mass-ratio threshold ^ = 0.1 (dotted lines). Mergers 
have limited effect at z > 2, but hy z = 0.1 they have 
dramatically boosted the space density of black holes with 
Mbh > 10^ Mq, to a level clearly inconsistent with the SWM 
and Vika et al. (2009) estimates. The discrepancy would be 
more severe if we did not include mass-dependent e in the 
G{z, Mbh) model (see Fig. llOp. We caution, however, that 
in this regime the local mass function estimates rely largely 
on extrapolation of the Mbh — cr* or Mbh — M« correlations 
into a range with limited observational constraints. This is 
also a mass range where typical galaxy hosts are gas poor, 
and it is not clear that "dry" galaxy mergers necessarily lead 
to black hole mergers. 

In general, the inclusion of mergers makes it even harder 
to fit the observed steep decline of black hole abundance at 
high masses. 



5.4 A A-dependent bolometric correction 

The inclusion of mergers puts otherwise acceptable models 
at risk of overpredicting the high mass end of the local black 
hole mass function. We now discuss one possible resolution 
of this tension, a A-dependent bolometric correction. 
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Figure 14. Binned data on the 2 — 10 keV X-ray bolometric 
correction as a function of Eddington ratio from Vasudevan & 
Fabian (2007) . The solid line is the analytical approximation used 
in this paper, while the long-dashed and dot-dashed lines are the 
results from Lusso et al. (2010). 

Following previous work (e.g., Elvis et al. 1994; Mar- 
coni et al. 2004; Shankar et al. 2004; Hopkins et al. 2007), 
we have so far assumed that the bolometric correction de- 
pends only on bolometric luminosity (all the X-ray data in 
SWM were converted to bolometric luminosities using the 
L-dependent bolometric correction by Marconi et al. 2004). 
However, some studies (e.g., Dai et al. 2004; Saez et al. 2008; 
and references therein) show signs for variations in the spec- 
tral energy distributions of AGNs with either redshift or ac- 
cretion properties. In particular, Vasudevan & Fabian (2007, 
2009) suggest that variations in the disc emission in the 
ultraviolet may be important to build the optical-to-X-ray 
spectral energy distributions of AGNs. From a sample of 
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Figure 12. Similar format to Figures [5] and [B] Comparison among the G(2,Mbh)i G{z) and G+P(2,Mbh) predicted Eddington ratio 
distributions {solid, red long-dashed, and cyan dashed lines, respectively) compared to the K06 data at z = 2 and 2 = 1, respectively 
{top and middle panels), and to H09 at z = 0.5 {lower left) and KH09 in the local Universe {lower right). In the lower panels, we also 
include curves for the steep power-law model introduced in Section 14.41 



54 AGNs from the Far Ultraviolet Spectroscopic Explorer 
(FUSE) and X-ray data from the literature, they claim evi- 
dence for a large spread in the bolometric corrections, with 
no simple dependence on luminosity being evident. Their 
results suggest instead a more well-defined relationship be- 
tween the bolometric correction and Eddington ratio, with 
a transitional region at an Eddington ratio of ~ 0.1, below 
which the X-ray bolometric correction is typically 15-25, and 
above which it is typically 40-70. As shown in Fig. 1141 we 



approximate their results by setting 

r 18 if A sC 0.105 

A'2-iokov(A) = <^ 54.85 + 26.78 logLx- 

[ -11.11 (log Lx)^ if 0.105 s; A 1 

(22) 

Our analytic approximation is shown as a solid line in 
Fig. 1141 against the binned data by Vasudevan and Fabian 
(2007; the fit and the data extend to A > 1), while the dot- 
dashed and long-dashed lines are more recent fits from Lusso 
et al. (2010) derived from a larger sample. As the latter au- 
thors point out, while a trend of bolometric correction with 
A might exist, determining its exact slope is still challenging 
given the large dispersion in the data. In this section we will 
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Figure 13. Predicted black hole mass function for the G(2, Mbh) model without mergers (dashed lines), with only major mergers [solid 
lines), and with minor and major mergers {dotted lines). The set of lines, from bottom to top, are the model predictions at 2: = 4, 2, 0.1, 
respectively. The data are as in Fig. |3] The cumulative effect of mergers is minor at z > 2 but becomes significant at low redshifts and 
high masses. 
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Figure 15. Influence of a A-dependent bolometric correction on 
the predicted black hole mass function. The rfot-dashed line shows 
the 2 = mass function of the G(z, Mbh) model with our stan- 
dard bolometric correction. The solid lines show the evolving 
mass function at 2 = 2, 1, and for the same model assuming 
the A-dependent bolometric correction of Eq. 1221 The dotted line 
is the A-dependent G(2, Mbh) model inclusive of mergers (with 
^ > 0.3). Data points and grey band are the same as those in 
Figs. [3l and [TOl 



use Eq.[22]as a reference, noting that the Lusso at al. (2010) 
fits provide consistent resuhs in the range of interest here. 

In our numerical formalism described in Appendix lAl 
it is straightforward to insert a A-dependent bolometric cor- 
rection. Given that the relation L oc Ar2-iokcv(A) x Lx oc 
AMbh, it implies that Lx oc [A/i5r2-iokcv(A)]AfBH. There- 
fore, having a A-dependent bolometric correction is equiv- 
alent to running the code replacing bolometric luminosi- 
ties with Lx- We thus solve Eq. [8] for computing the 
duty cycle by replacing the bolometric luminosity function 
^{L,z) on the left-hand side with the X-ray luminosity 



function ^x{Lx , z), and using effective Eddington ratios 
A' = A/A'2-iokcv(A). 

Most predictions of our models are not sensitive to 
this change of bolometric correction. In particular, we have 
checked that after taking into account the A-dependent con- 
versions between sample flux limits and bolometric lumi- 
nosities, both duty cycles and Eddington ratio distributions 
showed similar behaviours to the ones predicted by the 
G(z,Mbh) model with a luminosity-dependent bolometric 
correction. 

Nevertheless, the A-dependent correction does have a 
significant effect on the low- 2; black hole mass function. The 
soUd lines (at 2: = 0, 1, 2, from top to bottom) in Fig. llSl show 
the G{z, Mbh) model with a A-dependent bolometric correc- 
tion as in Eq. 1221 High-mass black holes have preferentially 
lower A in this model, and the lower bolometric luminosity at 
a given X-ray luminosity reduces the inferred growth rate of 
these massive black holes. Since G{z, Mbh) previously agreed 
well with the z = mass function, it now underpredicts the 
high mass end. 

A A-dependent bolometric correction would reduce 
the need for mass-dependent e in our G(2,AfBH) and 
G+P(2,Mbh) models (Section [STU Fig. [TO}, which was in- 
ferred partly from the high mass end of the local mass func- 
tion, though it was also supported by direct empirical ev- 
idence for mass-dependent e. Alternatively, a A-dependent 
bolometric correction could compensate the impact of merg- 
ers on our standard G{z, Mbh) model (see Fig. I13p to yield 
improved agreement with local mass function estimates. The 
dotted line in Fig. [15] is the predicted 2 = black hole 
mass function for the G(2,Mbh) model including a mass- 
dependent e, A-dependent bolometric correction, and black 
hole mergers with ^ > 0.3. Agreement with observational 
estimates is significantly improved relative to the standard 
bolometric correction case (uppermost solid line in Fig. llOp . 
Uncertainties in bolometric corrections — their normaliza- 
tion and their dependence on luminosity, A, or other factors 
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— remain an important source of uncertainty when testing 
evolutionary models of the black hole population against the 
local census of black holes. 



5.5 Specific Black Hole Accretion Rate 

Several authors have noted that the average black hole accre- 
tion rate has a redshift dependence morphologically similar 
to the cosmological SFR (e.g., Marconi et al. 2004; Merloni 
et al. 2004; Silverman et al. 2008a; Zheng et al. 2009; SWM). 
Fig. [16] adds a new piece of information to the co-evolution 
of black holes and galaxies, plotting the mean specific accre- 
tion rate of black holes of a given mass at a given redshift 
defined as 

(M(Mbh,^)) _ (A(AfBH,z));7(MBH,z) 

While the mean accretion rate does not depend on the -P(A) 
distribution but only on the assumed radiative efficiency, 
the mean specific accretion rate depends significantly on the 
input P(A), so it provides diagnostic power beyond that in 
the global rates. 

In the G{z) model (left panels), the specific mean ac- 
cretion rate declines with increasing black hole mass at all 
epochs. Cyan lines in the upper panels show the mean spe- 
cific star formation rate (SSFR) recently derived by Karim 
et al. (2011; see also Noeske et al. 2009) at the same redshifts 
(from their Table 4), with stellar masses simply re-scaled by 
a factor of 10~^ to convert to black hole masses. In this 
simple comparison with the G{z) model, the mean specific 
black hole accretion rate decreases with mass and increases 
with redshift in a remarkably similar way as the mean star 
formation rate, with a slope --^ A/g^'*. However, when we 
consider the 0(2, Mbh) model (right panels), which better 
matches data on black hole duty cycles and A-distributions, 
we predict higher specific accretion rates at high masses, re- 
flecting the higher duty cycles in this model. For simplicity, 
here we use a constant radiative efficiency, thus a constant 
ts in Eq. 1231 Silverman et al. (2009) also find that in rel- 
atively massive high-z galaxies the ratio between average 
black hole accretion rate and SFR is higher by up to an 
order of magnitude with respect to the classical 10""^, in 
fair agreement with the prediction of the G(z, Mbw) model, 
though dependent on the actual variations of radiative ef- 
ficiency with mass and/or time. In fact, we checked that 
including a mass-dependent radiative efficiency as in Eq. [T^] 
would line up the specific black hole accretion rate with the 
SSFR at all masses. 

The lower panels of Fig.[T6]show the mean specific black 
hole accretion rate as a function of redshift for different bins 
of black hole mass, as labelled. The grey bands indicate the 
uncertainties around the measured SSFR as catalogued and 
derived by Gonzalez et al. (2010) for galaxies with stellar 
mass ~ 5 X 10^ — 10^° Af©. The latter should be compared 
with only the specific accretion rate onto black holes with 
current mass ~ lO^Af0 (solid lines), but for completeness 
we also show the accretion rate for more massive black holes 
A^BH = IQ^Mq and A^bh = l(f Mq (long-dashed and dotted 
lines, respectively), as labelled. Overall, consistently with 
what is found in the upper panels, the models predict an 
accretion rate that tracks the SSFR, though the G{z) model 



tends to produce a specific accretion rate that steadily in- 
creases even at z > 2, at variance with the data for the star 
formation rate. The G(z,AfBH) model predicts a redshift 
dependence of the accretion rate morphologically similar to 
that of galaxies, a fact that might play some role in ex- 
plaining the still puzzling plateau at z > 2 of the galactic 
SSFR(2), which is poorly reproduced by semi-analytic mod- 
els and might require some extra source of early feedback 
(e.g., Weinmann et al. 2011). 



6 DISCUSSION 

Although we have considered a wide variety of models, with 
different P(A) shapes and different redshift and mass depen- 
dences of Ac and e, every one of these models shows signif- 
icant (factor of several) disagreement with at least one of 
the observational tests we have examined. This failure could 
indicate that our model assumptions are still too restrictive 
to describe the real black hole population — for example, 
we generally assume that the shape of P(A) is independent 
of redshift and black hole mass except for overall shifts in 
Ac, and we have considered restricted functional forms for 
these Ac trends and for P(A) itself. Alternatively, the prob- 
lem could lie in one or more of the data sets themselves, since 
these are frequently derived from noisy or uncertain estima- 
tors (e.g., for black hole masses) and from input samples that 
are subject to selection biases and incompleteness (e.g., for 
active galaxy fractions). Here we highlight aspects of the 
data that appear especially difficult to reproduce within our 
class of models, or where different data sets appear to drive 
the models in contradictory directions. Recall that all of our 
models reproduce SWM's estimate of the bolometric lumi- 
nosity function (summarized in Section 12. 2[) by construc- 
tion. SWM discuss remaining uncertainties in this luminos- 
ity function and their impact on inferred model parameters. 

The first tension within the data is at the low-mass end 
of the local black hole mass function (see Fig. [3]), where 
SWM (and a number of other studies) find a $(A'/bh) ris- 
ing to low masses but Vika et al. (2009) and some other 
studies (e.g., Graham et al. 2007) find a falling $(AfBH). 
This region of the mass function remains difficult to probe 
because of uncertainties in bulge-disk decomposition and be- 
cause the black hole mass correlations for spiral bulges are 
more uncertain than those for high mass, bulge-dominated 
galaxies. In our models, it is difficult to produce a turnover 
like that of Vika et al. (2009), though the G{z) model goes 
in this direction (Fig. llOp because it ascribes much of the 
low luminosity AGN activity to low-A accretion by massive 
black holes, hence reducing the growth of low mass black 
holes. However, the predicted mass function in this regime 
depends on the input luminosity function in a range that 
is largely extrapolated from brighter magnitudes. Therefore 
if "I>(A4"bh) really does turn over at low masses, a plausi- 
ble explanation is that the faint end of the AGN luminosity 
function is flatter than the one adopted here. 

A second tension, with more serious implications for 
the issues at the core of this paper, is between the nar- 
row Eddington ratio distributions measured by K06 and the 
broader distribution, peaking at lower Ac, measured by H09. 
The two data sets overlap in redshift, though H09 are at the 
low redshift end of the K06 range. H09, and KH09 at still 
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Figure 16. Mean specific black hole accretion rate (Eq. I23II predicted by the G(2:) (Ze/t) and G(2:, Mbh) {right) models as a function 
of black hole mass (top) and of redshift (bottom). In the top panels, cyan lines are the mean specific star formation rates as calibrated 
by Karim et al. (2011; their Table 4), with stellar masses simply scaled to black hole masses assuming a proportionality factor of 10""^, 
as measured in the local Universe. In the bottom panels, the grey area is the specific star formation rate with error bars calibrated by 
Gonzalez et al. (2010) for galaxies with stellar mass ~ 5 X 10^ — 10^^ Mq, which can be compared to the predictions for lO'^M© black 
holes (solid curves). 



lower z, favor a P(A) that is broad in shape (like G+P) and 
evolving to low Ac at low redshift, but it is difBcult to recon- 
cile such a model with the K06 histograms. The luminosity 
thresholds for the H09 and K06 data sets are very different, 
roughly Lmin = 10'*^ergs~'^ and Lmin = 10'*^erg s"'^, respec- 
tively, so there is no direct contradiction between the mea- 
surements. Possibly a model that allows a different shape 
of -P(A) for high and low mass black holes, or that allows 
a power-law that is steeper or further offset from the log- 
normal peak, could be made consistent with both sets of 
observations. The two samples might differ for more pro- 
found physical reasons, as recently emphasized by Trump et 
al. (2011), who showed that broad lines might disappear at 
A < 0.01 because of a change in accretion flow structure. 

The observed P(A) distributions also include errors in 
black hole mass estimations, so the intrinsic distributions 
should in principle be even narrower. Correcting for this ob- 
servational broadening would exacerbate the tension with 
our G{z,Mbh) and G-|-P(2,Mbh) models, which already 
predict broader P(A) distributions than those found by K06. 

Another important tension arises between the low duty 
cycles found for low redshift quasars by Schulze & Wisotzki 
(2010) and the much higher active galaxy fractions found by 



Kauffmann et al. (2003), with a gap that is roughly two or- 
ders of magnitude. The difference between Type 1 and Type 
2 AGN could possibly explain a factor of 10 — 20 (Greene 
& Ho 2009), though this factor is already large compared 
to conventional estimates of obscured-to-unobscured AGN 
ratios, and other empirical and physical effects might need 
to be invoked to explain such a strong discrepancy (e.g.. 
Trump et al. 2011). By adopting the Greene & Ho (2009) 
ratios in our predictions, we find models that are roughly 
consistent with both Kauffmann et al. (2003) and Schulze & 
Wisotzki (2010), but even our best cases disagree with one 
of these data sets by a factor of several in some black hole 
mass range (see Fig. [7| . We have converted Kauffmann et 
al.'s [Oin] luminosities to bolometric luminosities assuming 
a constant bolometric correction, and scatter or biases in this 
correction (see, e.g., Capetti 2011) might account for some 
of the discrepancy with Schulze & Wisotzki (2010). Best et 
al. (2005, see their Fig. 2) also find high AGN fractions for 
SDSS galaxies — 20-40% for L[oiii] > IO^'^-Lq (logL > 42) 
— comparable to those of Goulding et al. (2010). 

A fourth tension arises from the trend of higher ac- 
tive fractions for more massive galaxies found by Xue et al. 
(2010) and Bundy et al. (2010). If P(A) does not evolve, then 
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matching observed AGN luminosity evolution leads to down- 
sizing, i.e., a decrease of duty cycle with increasing black 
hole mass (Fig. |3]). A \c{z) that declines towards low red- 
shift can soften this trend, but within our considered range 
of models, it does not eliminate downsizing entirely. We have 
been able to produce a trend of rising duty cycle with ris- 
ing black hole mass by making Ac decrease towards higher 
Mbh, but the resulting models then tend to be inconsistent 
with the K06 Eddington ratio distributions, which do not 
show such a trend fFig. I12|) . Furthermore, the model trends 
of (7(Mbh) remain flatter than those found by Xue et al. 
(2010) and Bundy et al. (2008), and these models still yield 
falling U{Mbw) a.t z — 0.5. One caveat is that we are trans- 
lating the observed galaxy stellar masses to corresponding 
black hole masses assuming a linear relation with no scatter. 
If the scatter between black hole mass and galaxy mass were 
large at these redshifts, then the trend between active frac- 
tion and galaxy mass could be partly induced by a higher 
probability for black holes of fixed mass to be active if they 
reside in more massive galaxies. 

In Section 14.41 we considered a model in which P{\) 
rises steeply towards low A (with P oc A""'^), which helps 
produce a rising ?7(Mbh) trend for a luminosity thresholded 
sample because more massive black holes can radiate at low 
A while remaining above threshold (Aird et al. 2011; Mainieri 
et al. 2011). We restricted this steep P(A) to z < 0.7, since 
at higher redshifts it leads to duty cycles above unity. With 
this model (inspired by that of Aird et al. 2011), we are 
able to obtain a rising U{AiBH) at z = 0.5. However, the 
prediction for the observed P{X) disagrees with the data of 
K06, H09, and KII09, a discrepancy also noted by Aird et 
al. (2011). 

A final tension arises at the high end of the black hole 
mass function, where the models tend to overpredict the 
data. While models with a mass-independent P(A) are ac- 
ceptable within the estimated observational uncertainties, 
adding mass dependence of Ac to produce a rising U{Mbu) 
leads to a substantial overprediction of the abundance of 
the most massive black holes. We have mitigated this prob- 
lem in our G{z, Mbh) model by also increasing the radiative 
efficiency at high Mbh (Fig. llOf) . This solution has some ob- 
servational support, as discussed in Section [431 On the other 
hand, including mergers as calculated in Section [5 . 31 worsens 
the overprediction of $(Mbh) at high masses (Fig. I13|) . 

There are other factors that may ameliorate the dis- 
crepancy between models and observational determinations 
of the abundance of the most massive black holes. A A- 
dependent (or simply lower) bolometric correction may re- 
solve this tension, as discussed in Section 15.41 The model 
shown by the dotted curve in Fig. 1151 which includes merg- 
ers, is acceptable within current uncertainties. The high end 
of the mass function relies on extrapolation of the Mbh- 
bulge relations into a regime where there are few calibrat- 
ing galaxies (Schulze & Gebhardt 2011). Another effect that 
can help fitting the AGN luminosity function with a steeper 
decline of the black hole mass function at high masses is 
anisotropic emission by AGNs, which would give rise to ap- 
parent values of A (i.e., inferred from the observed flux in 
our direction under the assumption of isotropy) that are 
occasionally larger than unity. The majority of the most lu- 
minous AGN would then correspond to objects that have 
their brightest direction of emission pointing to us, rather 



than to AGN with the highest black hole masses, and the 
mass accretion rates of these AGN would be lower than their 
observed luminosities suggest. (Note that in our models we 
have imposed a cutoff on the P(A) distribution at A > 1.) 

If, in light of this discussion, we take a rather gen- 
erous view of the systematic uncertainties in the observa- 
tional constraints we have considered, then our G{z, Mbh) 
and G-fP(2,MBH) reference models may be viewed as at 
least moderately successful, while the other reference mod- 
els — S, G, G-l-P, and G(z) — all fail drastically on at 
least one observable. Table [2] summarizes the observational 
comparison based on an admittedly subjective assessment 
of the resuhs shown in Figures El 13 El H [101 and [H We 
consider as constraints the SWM estimate of the local black 
hole mass function, the high-z Eddington ratio distributions 
from K06, the low-z Eddington ratio distributions from H09 
and KH09, and the duty cycle estimates from active galaxy 
fractions at 2; « 0.2 shown in Figure [71 Recall that all models 
reproduce our input AGN luminosity function by construc- 
tion. We assign a / when a model reasonably describes an 
observation with some allowance for systematic uncertainty, 
an X when it clearly fails, and a — for intermediate cases. 
The non-evolving models and G(z) model all fail to match 
the low-2 -P(A) or the low-z ?7(Mbh, imin), or both. The 
G(2,A/bh) and G+P{z, AIbu) models have no such dras- 
tic failures, though they are far from perfect matches to 
the data. However, none of our models reproduce the trend 
of duty cycle with black hole mass illustrated in Figure [HI 
though the G(z, Mbh) model is the least discrepant. 



7 CONCLUSIONS 

We have extended the formalism of continuity-equation 
modeling of the black hole and AGN populations to allow 
a distribution of Eddington ratios -P(A). With this broader 
class of models we have addressed two new categories of ob- 
servations, direct estimates of Eddington ratio distributions 
of active black holes and estimates of duty cycles from active 
galaxy fractions. Both of these categories have been areas of 
intense observational investigation over the last five years, 
and as representative examples we have concentrated on 
P(A) estimates from KoUmeier et al. (2006; K06) at 2 ^ 1, 
Hickox et al. (2009; H09) a.t z x 0.5, and Kauffmann & 
Heckman (2009; KH09) at 2 « 1, and on active galaxy frac- 
tion data from Bundy et al. (2008) and Xue et al. (2010) 
at z ^ 0.5 and from Goulding et al. (2009), Kauffmann et 
al. (2003), and Schulze & Wisotzki (2010) at low redshift. 
We account for the effective luminosity thresholds of these 
analyses in our model predictions, though the thresholds are 
not always clearly defined. 

As forms for P(A) we consider a Gaussian in log A (G) 
and a Gaussian with a power-law extension to low- A (G-l-P). 
If the radiative efficiency e and characteristic Eddington ra- 
tio Ac (where the Gaussian peaks) are held fixed, then chang- 
ing from a single A value to either of these distributions has 
little impact on the evolution of the black hole mass function 
inferred from continuity-equation modeling. As in the single- 
A models of SWM (and similar models by Shankar et al. 2004 
and Marconi et al. 2004), we find that models incorporat- 
ing the observed AGN luminosity function and parameter 
values e ~ 0.07 and Ac ~ 0.25 yield a good match to ob- 
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Model 



'S>bh{Mbh,z), z = P{X\Mbh,z), z>0.5 P{X\Mbu, z), z ^ 0.5 U (Mbh , L^in , z) 



S 
G 

G+P 
G(z) 

G(z,Mbh) 
G+P(z,Mbh) 



/ 
/ 

/ 
/ 



X 
/ 
/ 
/ 

/ 



X 
X 
X 
X 

/ 



X 
X 
X 

/ 
/ 



Table 2. List of the Reference Models listed in Table [T] along with a qualitative assessment of their agreement with the data. We assign a 
/ when a model reasonably describes an observation with some allowance for systematic uncertainty, an X when it clearly fails, and a — 
for intermediate cases. The P(A|AfBHi z) at 2 > 0.5 column refers to the K06 data (Figures [5l and [T2I 1. The P(A|Mbhi 2) at z ^ 0.5 refers 
instead to the H09 and KH09 data (Figures [6] and [12} . Finally, the (7(Mbh , ^mim z) column refers to the multiple data sets reported in 
Figures [3 and [8] 



servational estimates of the black hole mass function in the 
local Universe. The predicted black hole duty cycle declines 
rapidly with decreasing redshift at z < 2, and at redshifts 
2 ^ 1 it declines sharply with increasing black hole mass 
over the range lO'^M© — W^Mq ("downsizing" evolution). 
An unattractive feature of the G+P model is that it requires 
a high space density of massive black holes already present 
at z = 6, since only the small fraction of black holes with 
high A are luminous enough to contribute to the observed 
luminosity range at this redshift. This large "seed" popula- 
tion appears physically unrealistic, and the low implied duty 
cycle for high- luminosity quasars contradicts evidence from 
their strong observed clustering at 2; ~ 4 (Shen et al. 2007; 
White et al. 2008; Shankar et al. 2010b). We conclude (in 
agreement with Cao 2010) that any extended low-A tail of 
P(A) must develop at lower redshifts (z 3) rather than 
being a redshift-independent feature of black hole fueling 
processes. 

Our models with redshift-independent P(A) predict 
duty cycles at 2 ^ 1 that are far below observational es- 
timates from active galaxy fractions, and they do not match 
the \ow-z Eddington ratio distributions estimated by H09 
and KH09. Motivated by these discrepancies, we introduce 
the G{z) model with a redshift-dependent Ac (Eg.llSp. shift- 
ing P{X) to lower Eddington ratios at low redshift. This 
model yields better agreement with the observations. The 
low-redshift P{X) remains narrow compared to the H09 and 
KH09 estimates, and the model still predicts a duty cycle 
that declines with increasing black hole mass, in contra- 
diction to the Bundy et al. (2008) and Xue et al. (2010) 
finding of higher AGN fractions in more massive galaxies. 
We therefore introduce an additional mass dependence of 
Ac (Eq. 1161 Ac oc Afgn '^) at each redshift, which, relative 
to mass-independent models, maps massive black holes to 
lower luminosity, more numerous AGN. At the same time, 
we introduce a steady broadening of the Gaussian P(A) to- 
wards low redshift (Eq. ll7p . or, for the G+P model, a steady 
drop in the minimum Eddington ratio (Eg. llSfl that removes 
the need for an unrealistic seed population at z = 6. 

The predictions of the resulting models, G(z, Mbh) and 
G+P(2:, A/bh), are fairly similar. Both models achieve a rea- 
sonable match to the H09 and KH09 Eddington ratio dis- 
tributions, though they underpredict the KH09 distribution 
at A ^ 10"'^'^. They produce approximate agreement with 
the low-redshift duty cycle estimates if we adopt the large 
(factors of 10 - 20) ratios of Type II to Type I AGN advo- 
cated by Greene & Ho (2009), though there are still factor 



of several discrepancies for some data sets at some A^bh val- 
ues. Both models predict duty cycles that are flat or weakly 
rising with black hole mass at 2: > 1, thus improving the 
agreement with Bundy et al. (2008) and Xue et al. (2010), 
though neither model can reproduce the Xue et al. (2010) 
trend at 2 = 0.5. 

Our models with mass-dependent P(A|AfBH,2) exhibit 
tension with the K06 Eddington ratio distributions, espe- 
cially at 2 « 1, since the redshift and mass dependence of Ac 
drive the predicted -P(A) distributions to peak at log A ~ — 1 
while the observed distributions peak at log A ~ —0.6. 
Adding a steep (P cx A"'''^) power-law at late times, simi- 
lar to the model of Aird et al. (2011), can better match the 
Xue et al. (2010) duty cycle data ai z — 0.5, but it also 
spoils the match with the H09 and KH09 Eddington ratio 
distributions. 

Boosting the duty cycle, and thus the growth, of mas- 
sive black holes tends to overproduce the high-mass end of 
the local mass function relative to observational estimates. 
Our standard versions of the G(2, A/bh) and G-|-P(2, A4"bh) 
models therefore incorporate a mass-dependent radiative ef- 
ficiency (Eq. [191), following Cao & Li (2008). The higher 
efficiency assumed for higher mass black holes reduces their 
inferred growth, restoring agreement with the local mass 
function. Alternatively, a A-dependent bolometric correction 
(Vasudevan & Fabian 2007) can lower the inferred growth 
of massive (hence lower A) black holes in these models, ob- 
viating the need for mass- dependent e. However, black hole 
mergers at the rate suggested by black hole merger statis- 
tics also raise the high mass end of the mass function at 
low redshifts, so in a complete model the mass-dependent 
e may be be required even with the A-dependent bolomet- 
ric correction. Furthermore, Davis & Laor (2011; but see 
also Raimundo et al. 2011 and Laor & Davis 2011) have 
presented direct evidence for mass-dependent e from quasar 
spectral energy distributions, and such a dependence could 
also explain the discrepancy between the e > 0.2 inferred 
for luminous, strongly clustered quasars at z — 4 (Shankar 
et al. 2010b) and the average e ~ 0.07 implied by matching 
the local black hole mass density (SWM and numerous refer- 
ences therein). Each of these arguments for mass-dependent 
e rests on somewhat shaky ground, but they all point in 
the same direction. At 2 > 1, the black hole mass functions 
implied by our models are fairly insensitive to mergers and 
only moderately sensitive to other model assumptions. 

In agreement with previous studies, we find that the 
growth of the black hole population tracks the overall evolu- 
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tion of the cosmic star formation rate. We also find generally 
good agreement between the mean specific black hole growth 
rates (Mbh)/Mbh and Karim et al.'s (2011) estimates of the 
specific star-formation rates of star-forming galaxies over the 
full range 0.3 ^ z ^ 3, if we simply translate stellar masses 
to black hole masses with a constant scaling factor of 10""^. 

Over the last five years, measurements of AGN cluster- 
ing have improved dramatically in precision, redshift extent, 
and luminosity range. These measurements provide valuable 
constraints on the relation between active black holes and 
their host dark matter halos, which we have previously ex- 
plored in the context of single-A accretion models (Shankar 
et al. 2010b, c). We will extend our clustering studies to in- 
clude P(A) distributions and mergers in future work. 

With P{X) and radiative efficiency allowed to depend on 
redshift and black hole mass, our models have become rather 
elaborate despite their simple physical basis. This complex- 
ity reflects the growing richness of the observational data, es- 
pecially the measurements of Eddington ratio distributions 
and active galaxy fractions over a wide range of redshift, 
mass, and luminosity, from a variety of data sets. In order 
of decreasing robustness, our key qualitative conclusions are 
(a) that the characteristic Eddington ratio Ac declines at 
low redshift, (b) that the P(\) distribution broadens at low 
redshift, (c) that more massive black holes have lower Ac, 
and (d) that more massive black holes have higher radiative 
efficiency. 

Despite the flexibility of our models, and despite inves- 
tigating many variants beyond those discussed in the paper, 
we have not found a model that fully reproduces all of the 
observational constraints we have considered. The remain- 
ing discrepancies presumably reflect some combination of 
inadequate models and systematic errors in the data sets, 
as discussed in detail in Section [6] On the model side, we 
have assumed restricted functional forms for the mass and 
redshift dependence of -P(A), and more general behavior or 
sharper evolutionary transitions may be required to match 
the data. Observationally, estimates of black hole masses and 
bolometric luminosities are both subject to systematic un- 
certainties, and even when these estimates are correct in the 
mean, scatter can have important effects on inferred trends 
and distributions. Many of the tensions between the mod- 
els and the data and among the data sets themselves revolve 
around the seemingly disparate trends found for optically lu- 
minous, broad-line quasars and varieties of low luminosity. 
Type II AGN. We have followed standard practice in re- 
lating these two populations by a simple obscuration factor, 
but a more nuanced relation between the different categories 
of active black holes (e.g.. Trump et al. 2011) may be crucial 
to resolving some of the tensions highlighted here. 

Continuity-equation models draw on the inevitable link 
between luminosity and mass accretion to tie the observable 
population of AGN to the evolving population of supermas- 
sive black holes that power them. They provide a powerful 
framework for linking empirical studies that probe a vari- 
ety of observables across a wide span of redshift, luminosity, 
and black hole mass. As these empirical studies continue to 
improve in precision, dynamic range, and control of system- 
atic uncertainties, they will refine the models into a tightly 
constrained history of the cosmic black hole population. 
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APPENDIX A: SOLVING THE CONTINUITY 
EQUATION WITH BROAD INPUT 
EDDINGTON RATIO DISTRIBUTIONS 

There is not a unique way to solve Eq. [S] Steed & Wein- 
berg (2003) adopted an input parametrical double power-law 
duty cycle, and tuned the parameters to predict the AGN 
luminosity function. More recently, a similar technique has 
been adopted by Cao (2010), who performed a detailed 
to up-to-date AGN luminosity functions at all redshifts, and 
then solved the continuity equation having full information 
on the input duty cycle. 

Following these previous works, and in particular Cao 
(2010), at any redshift z we parameterize the input active 
mass function A^act(AfBH, z), i.e., the product of the duty 
cycle and the black hole mass function, by a double-power 
law of the type 

iVact(MBH, Z) = -J . (Al) 

At all redshifts, we then determine the parameters A^o, ce, P, 
and by first computing the convolution with the input 

Eddington ratio distribution P(A|A/bh,z) and then finding 
the best-fit to the AGN luminosity function 

$(L,2) = J d\ogXP{X\MBH,z)N^ct{MBu,z). (A2) 

Having A^act {Mbh , z) at all redshifts we then compute 
the average accretion rate from Eq.[S]and then compute the 
black hole mass function from the continuity equation in 
Eq. [1] The duty cycle at all redshifts and black hole masses 
is then simply computed from the ratio of A'act(MBH, z) and 
the total black hole mass function. The method described 
above, which we take as a reference throughout the paper, 
allows one to describe the accretion histories of black holes 
for any continuous input Eddington ratio distribution. How- 
ever, it also relies on assuming an a priori shape for the active 
mass function of black holes. 

We have developed another method to solve for the duty 
cycle in Eq. [8] that does not rely on any a priori shape for 
the active mass function. We first assume that the active 
black holes of mass within A/bh at a given redshift z with 
duty cycle ?7(Mbh, z) accrete following a predefined discrete 
Eddington ratio distribution Pj(MBH, z). We take A'' values 
of the Eddington ratio A = Aj, with j = 1, .., A*', and set the 
relative probability for a black hole to accrete at Aj equal to 
Pj(MBH,2). The fraction of active black holes accreting at 
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A = \j at redshift z is then given by Pj{MBH, z)U{AIbh, z), 
and the distribution Pj must be such to satisfy the condition 

J2 Pj{MBH,z)UiMBn, z) = (/(A/bh, z) (A3) 



so that the total fraction of active black holes of mass Mbh 
is given by )7(AfBH,2). We then solve Eq. ((HI numerically 
by computing the duty cycle at each point A/"bh in the grid 
as 

C/(Mbh,2) = {P,^,^^^BB{MBYi,z))-^ [$(L,2)|z,«A,.^^^^Mbh ^ 
- 5Z 'Pj*BH(AfBH,2:)|£o,AjM^jj;M^jj>J\/BH] 

A black hole of mass AfBH can radiate at the maximum lu- 
minosity L oc y^jMAx -^BH set by the maximum assumed Ed- 
dington ratio ^^jMAx ■ More massive black holes with mass 
Afgjj > A^BH can in principle radiate with similar lumi- 
nosity L but at lower Eddington ratios Aj < ^j^Ax i ac- 
cording to the relation L cx Aj A/bh- Therefore Eq. A4 
computes the duty cycle at a given bin around ATbh sub- 
tracting from the luminosity function computed at the bin 
L oc Ajj^^x-^BH, the contributions of more massive black 
holes with A^bh > AfsH shining within the same bin around 
L with lower Eddington ratios. 

Note that Eq. A4 is implicit, as in order to compute 
[/(A^bh, z) it requires knowledge of the duty cycle at higher 
masses U{Mbii > AfBH,^). However, the more luminous 
sources in the AGN luminosity function with luminosity 
Lmax can only be produced by the more massive black 
holes in the grid with mass A^bh.max accreting at the max- 



imum Eddington ratio \j 

?7(AfBH,MAX, z) 



Therefore we can set 



$(L,2) 



'&bh(A/bh,max, z) 



■^iMAX ^"BH.MAX 1 

(A5) 

and then iteratively solve Eq. A5 for lower and lower mass 
black holes. We have checked that our results are not depen- 
dent on the exact choices for the maximum black hole mass 
in the grid A^bh.max or the maximum luminosity Lmax, 
as long as they are sufficiently large (e.g., logA/BH,MAx > 
9.5 — 10 and logLAfAx > 48). For all our models we set to 
^Jmax ~ 1- 

We have checked that the discrete method described 
above yields consistent results with the continuous method 
when the same input P(A|A4"bh,2) distribution is adopted. 
However, we do not use the discrete method as a reference 
for several reasons. Being discrete, it generates oscillations 
in the predicted black hole mass function and duty cycle, 
absent from the continuous model. Moreover, due to its dis- 
crete nature, the former technique sometimes yields discon- 
tinuities at late times, especially when a complicated, mass- 
dependent, mapping between luminosity and halo mass is 
assumed. Other methods to solve Eq. [T] in the presence of 
broad Eddington ratio distributions have been discussed by, 
e.g., Cao & Li (2008), although those techniques, based on 
converting luminosity-based P{X) to mass-based ones, work 
well mainly for Gaussian distributions. 

In Section [5.41 we discuss models characterized by a A- 
dependent bolometric correction. For this class of models the 
solution to the continuity equation is still straightforward 
even in the presence of broad P(\\Mb¥i, z) distributions. 
The same exact procedure detailed above can be extended 



to A-dependent bolometric corrections by making use of the 
equality Lboi ~ L2-iokcvK{\) — AI/Edd and simply assum- 
ing a new scale of "effective" A' = X/K{X) that map X-ray 
luminosities L2-iokoV to Eddington luminosities, and using 
the X-ray luminosity function into Ea. lA2l instead of the full 
bolometric luminosity function. 



APPENDIX B: THE MEAN RELATION 
, BETWEEN BLACK HOLE AND HOST HALO 

■ Our models predict, by construction, the duty cycle and 
black hole mass function at all times during the evolution. 
It is therefore possible to infer the mean relation between 
black hole mass and halo mass via the cumulative relatiorjfl 
(e.g., Haiman & Hui 2001; Martini & Weinberg 2001; Wyithe 
& Loeb 2005) 



d log Af$i^ ( AT, z) 



dlogAfBH$BH(AfBH,2) (Bl) 

log Mbh 

with both ^BH and in units of comoving Mpc~^dex~^ 
for i?o = 70 km s" ^ Mpc" \ 

The halo mass function $h is "corrected" to include 
subhaloes as 

$HiM,z)=$ST{M,z)+ J nsiiiM,M')$sT{M',z)d\ogM' , 

(B2) 

where >I>sT(Af, z) is the Sheth & Tormen (1999j!| halo mass 
function and nsH(Af, Af') is the mass function from Giocoli, 
Tormen & van den Bosch (2008) which provides the number 
of subhaloes with unstripped mass within M and M + dM 
contained by haloes of mass M' and M' + dM' . Eq. B2 
significantly steepens the halo mass function below M ~ 
lO^^h-^M0,asit increases the number of lower mass haloes 
more, thus somewhat altering the mean correlation between 
haloes and black holes in this mass regime. 



APPENDIX C: COMPUTING BLACK HOLE 
MERGERS 

In the presence of mergers two additional terms should be 
added on the right hand side of Eq. [T] 

S^n = \ C di(^^{i,M)n^[M{M^^,z),z^ '^^ 



At 



At 



-{^,M)nh [Af (A/^'h,^),^] 



dAfBH 

dM 

dMBH 



° In this work we do not consider any scatter in the median black 
hole-halo relation. In future work we will probe the most suitable 
black hole-halo relation against AGN clustering, scatter, subhalo 
accretion histories. Changing the ATbh — A/ relation would alter 
the impact of mergers relative to that discussed in Section 15.31 
but we expect such changes to be small. 

^ More recent and detailed analysis of the halo mass functions 
have been performed (e.g.. Tinker et al. 2008); however, we here 
still use the Sheth & Tormen recipe to make contact with previous 
works and note that the exact choice of halo mass function does 
not alter our conclusions. 
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is the merger rate of incoming smaller mass black holes with 
mass Mbh = MbhC/(1 + and Mbh = Mbh/(1 + £,) that 
merge into a black hole of final mass Mbh, and 




is the merger rate of black holes with initial mass A/bh 
that merge into more massive black holes of mass A/bh = 
A/bh(1 + and A/^j'h = A/bh(1 + £,)■ In both Eqs. CI 
and C2 we set ^min = 0.1 or ^min = 0.3, and add the factor 
of 1/2 to avoid double counting; an additional factor of 1/2 
is present in Si-n to take into account that two black holes 
merge into one. 

The probabilities of halo mergers per unit time as a 
function of mass ratio, redshift, and remnant mass are taken 
from Fakhouri & Ma (2008). Following Shen (2009), we then 
insert a delay between host dark halo and host galaxy merger 
rate determined by the dynamical friction time of typical 
haloes of that mass. The probability of black hole mergers 
per unit time is then simply given by the "delayed" halo 
merger rate, i.e., 

Pu...^^>. M)nn[M{M^^, z),z]^ B^[M, g(z), z] , (C3) 

where z' is the actual epoch of the merger of the progenitor 
haloes (see Shen 2009 for details). By simply knowing, at 
each timestep, the mapping between infalling halo mass and 
its central black hole (given by Eq. [Bl| ). we can then com- 
pute the expected average rate for any black hole merger 
event. 



APPENDIX D: A COMPREHENSIVE 
COLLECTION OF ACTIVE GALAXIES AT 
DIFFERENT MASSES AND REDSHIFTS 

In this Appendix we present an overview of a detailed liter- 
ature search aimed at extracting an approximate estimate 
of the duty cycle of active galaxies at redshifts z > 1 above 
some given luminosity. 

At z ~ 2 several estimates of the number density 
of certain type of galaxies, along with their AGN frac- 
tion, have been measured. Caputi et al. (2006) estimated 
the cumulative number density of /(T-selected galaxies at 
high redshifts from the GOODS-Chandra Deep Fields to 
be < 2 X 10~*Mpc~^. They also found that the majority 
of the most massive galaxies in their sample, with stellar 
masses > 2.5 x 10^^ A/©, are ultraluminous infrared galax- 
ies, out of which at least 15% show AGN signatures from 
their X-ray luminosities. They claim the AGN fraction to 
decrease by about a factor of two moving to lower masses 
and lower infrared luminosities, although this may be partly 
due to selection effects (Caputi et al. 2006). 

In order to extract an absolute value of the AGN frac- 
tion from the AGN percentage revealed in the Caputi et 
al. (2006) sample, we need to correct for the fact that the 
density of galaxies in their sample may be below the ac- 
tual total number of galaxies of the same mass at the same 
redshift. Drory et al. (2005) derived the total galaxy stellar 



mass function at 2 < 2: < 2.5 from / and K selected galax- 
ies in the FORS Deep and GOODS Fields, respectively (see 
also Fontana et al. 2006). More recently, near-IR estimates 
come from Perez-Gonzalez et al. (2008) and Marchesini et 
(C2i)l. (2009). Given the large clustering of high redshift mas- 
sive galaxies (e.g., Foucaud et al. 2007; Magliocchetti et al. 
2007; Quadri et al. 2007) field-to-field variations may signif- 
icantly affect these results (e.g., Marchesini et al. 2007), but 
the nice agreement with the estimate by van Dokkum et al. 
(2006) derived from a large area is also reassuring. 

We thus find that the _K'-selected galaxy number den- 
sity derived by Caputi et al. (2006) amounts to about 60% 
of the cumulative number density of galaxies in the stellar 
mass function with mass above 10^^ Mq. We therefore take 
a mean of 0.15 x 0.6 ~ 0.1 as representative of the AGN 
fraction in galaxies with at 2: ~ 2 of this mass. 

Erb et al. (2006) find that only 5 out of 114 Lyman 
galaxies host AGNs, consistent with previous claims from 
similar samples by, e.g., Nandra et al. (2002, 2005), Gawiser 
et al. (2007). Adelberger & Steidel (2005) estimated an AGN 
duty cycle of about ~ 15% (although with large error bars) 
by combining clustering measurements and space densities 
of AGNs within their large sample of star forming galaxies 
at 1.8 < z < 3.5. Kriek et al. (2007) instead found from a 
sample of 20 -ft'-selected galaxies, that the fraction of AGNs 
inferred from several emission line diagnostics, was ~ 20% 
up to 40% for the more massive galaxies, and rapidly drop- 
ping at lower masses. The Kriek et al. results are at variance 
with the AGN fraction estimated in UV-selected galaxies by 
Erb et al. (2006) but in reasonable agreement with Caputi 
et al. (2006). However, Kriek et al. (2007) proved that their 
estimates are actually consistent with those from UV sam- 
ple. By binning the Erb et al. sample in stellar mass they 
found in fact that the AGN fraction in galaxies with stellar 
mass of 10^^ Mq is about 16% rising to ~ 40% at higher 
masses. Similar results were also claimed by Reddy et al. 
(2005), shown by a filled circle in the same Figure, who esti- 
mated that ~ 25% of their A's-selected galaxies in the deep 
Chandra fields are active. Note that the Kriek et al. and Erb 
et al. results are made somewhat uncertain by the fact that 
we do not have a proper estimate of the comoving abun- 
dance of galaxies in their samples. Papovich et al. (2005) 
also estimated an AGN fraction of 25% among a sample 
of distant red galaxies, which reduces to an overall 20% if 
distant red galaxies comprise about 70% of the overall pop- 
ulation of massive high redshift galaxies (van Dokkum et al. 
2006; Marchesini et al. 2007). 

The majority (> 75%) of the Submillimeter galaxies 
observed with deep SCUBA surveys at 850/i and ultradeep 
X-ray observations from the 2Msec Chandra Deep Field 
North, have been observed to host a central AGN (Alexan- 
der et al. 2005). The mean BH mass has been calibrated 
to be ~ 10* Mq and the mean Eddington ratio distribution 
peaked around 0.2 < A < 1 (Borys et al. 2005; Alexan- 
der et al. 2008). However, assuming the comoving density of 
SCUBA galaxies to be less than 1/4 the cumulative number 
density of Tf-selected massive galaxies (e.g., Chapman et al. 
2005), we get an average AGN fraction of ~ 17%, compara- 
ble to the result by Caputi et al. (2006). 

At lower redshifts estimates of the AGN fraction in large 
samples of galaxies have been recently estimated by Bundy 
et al. (2008), from DEEP2 and AEGIS surveys. In their sam- 



© 2012 RAS, MNRAS 000, [T]- ?? 



Accretion- Driven models of Black Holes 



pie AGNs have been detected in the 2-10 keV hard band of 
Chandra and therefore they are representative of the over- 
all AGN population (except for highly obscured, Compton 
thick AGNs). We find the Bundy et al. (2008) stellar mass 
function to overall in good agreement with other studies, 
within a factor of ~ 2. The ratio between the active and 
inactive stellar mass functions estimated by Bundy et al. 
(2008), yields a duty cycle increasing with black hole mass. 
The AGN fraction is ~ 10% raising to ~ 50% at z ~ 1.1, and 
from 0.01 to 0.1 at 2 ~ 0.5. This result is in good agreement 
with the estimates by Lehmer et al. (2007) who find that 
about 5% of the early-type galaxies in the Extended Chan- 
dra Deep Field-South contain AGNs in the local universe, 
increasing as (1 + z)^ at higher redshifts (i.e., about 17% 
at 2: ~ 0.5). A closer inspection by Silverman et al. (2008b) 
on an X-ray selected sample at redshifts 0.63 ^ z < 0.76, 
has confirmed an AGN fraction of (15± 5)%, decreasing to 
~ 2 — 6% for lower luminosity galaxies. Similarly, Shi et al. 
(2008) find a lower hmit of 2% for X-ray selected AGNs at 
2: ~ 0.5 in low mass galaxies increasing to ~ 10% in the more 
massive systems, in good agreement with what inferred by 
Bundy et al. (2008). The Bundy et al. (2008) average esti- 
mates also agree well with the Montero-Dorta et al. (2009) 
results on the AGN fraction of galaxies residing in different 
environments at 2: 1 from AEGIS and DEEP2. 
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